28class PseudoDifferentialOperator:
  29    """
  30    Pseudo-differential operator with dynamic symbol evaluation on spatial grids.
  31    Supports both 1D and 2D operators, and can be defined explicitly (symbol mode)
  32    or extracted automatically from symbolic equations (auto mode).
  33
  34    Parameters
  35    ----------
  36    expr : sympy expression
  37        Symbolic expression representing the pseudo-differential symbol.
  38    vars_x : list of sympy symbols
  39        Spatial variables (e.g., [x] for 1D, [x, y] for 2D).
  40    var_u : sympy function, optional
  41        Function u(x, t) used in auto mode to extract the operator symbol.
  42    mode : str, {'symbol', 'auto'}
  43        - 'symbol': directly uses expr as the operator symbol.
  44        - 'auto': computes the symbol automatically by applying expr to exp(i x ξ).
  45
  46    Attributes
  47    ----------
  48    dim : int
  49        Spatial dimension (1 or 2).
  50    fft, ifft : callable
  51        Fast Fourier transform and inverse (scipy.fft or scipy.fft2).
  52    p_func : callable
  53        Evaluated symbol function ready for numerical use.
  54
  55    Notes
  56    -----
  57    - In 'symbol' mode, `expr` should be expressed in terms of spatial variables and frequency variables (ξ, η).
  58    - In 'auto' mode, the symbol is derived by applying the differential expression to a complex exponential.
  59    - Frequency variables are internally named 'xi' and 'eta' for consistency.
  60    - Uses numpy for numerical evaluation and scipy.fft for FFT operations.
  61
  62    Examples
  63    --------
  64    >>> # Example 1: 1D Laplacian operator (symbol mode)
  65    >>> from sympy import symbols
  66    >>> x, xi = symbols('x xi', real=True)
  67    >>> op = PseudoDifferentialOperator(expr=xi**2, vars_x=[x], mode='symbol')
  68
  69    >>> # Example 2: 1D transport operator (auto mode)
  70    >>> from sympy import Function
  71    >>> u = Function('u')
  72    >>> expr = u(x).diff(x)
  73    >>> op = PseudoDifferentialOperator(expr=expr, vars_x=[x], var_u=u(x), mode='auto')
  74    """
  75
  76    def __init__(self, expr, vars_x, var_u=None, mode='symbol'):
  77        self.dim = len(vars_x)
  78        self.mode = mode
  79        self.symbol_cached = None
  80        self.expr = expr
  81        self.vars_x = vars_x
  82
  83        if self.dim == 1:
  84            x, = vars_x
  85            xi_internal = symbols('xi', real=True)
  86            expr = expr.subs(symbols('xi', real=True), xi_internal)
  87            self.fft = partial(fft, workers=FFT_WORKERS)
  88            self.ifft = partial(ifft, workers=FFT_WORKERS)
  89
  90            if mode == 'symbol':
  91                self.p_func = lambdify((x, xi_internal), expr, 'numpy')
  92                self.symbol = expr
  93            elif mode == 'auto':
  94                if var_u is None:
  95                    raise ValueError("var_u must be provided in mode='auto'")
  96                exp_i = exp(I * x * xi_internal)
  97                P_ei = expr.subs(var_u, exp_i)
  98                symbol = simplify(P_ei / exp_i)
  99                symbol = expand(symbol)
 100                self.symbol = symbol
 101                self.p_func = lambdify((x, xi_internal), symbol, 'numpy')
 102            else:
 103                raise ValueError("mode must be 'auto' or 'symbol'")
 104
 105        elif self.dim == 2:
 106            x, y = vars_x
 107            xi_internal, eta_internal = symbols('xi eta', real=True)
 108            expr = expr.subs(symbols('xi', real=True), xi_internal)
 109            expr = expr.subs(symbols('eta', real=True), eta_internal)
 110            self.fft = partial(fft2, workers=FFT_WORKERS)
 111            self.ifft = partial(ifft2, workers=FFT_WORKERS)
 112
 113            if mode == 'symbol':
 114                self.symbol = expr
 115                self.p_func = lambdify((x, y, xi_internal, eta_internal), expr, 'numpy')
 116            elif mode == 'auto':
 117                if var_u is None:
 118                    raise ValueError("var_u must be provided in mode='auto'")
 119                exp_i = exp(I * (x * xi_internal + y * eta_internal))
 120                P_ei = expr.subs(var_u, exp_i)
 121                symbol = simplify(P_ei / exp_i)
 122                symbol = expand(symbol)
 123                self.symbol = symbol
 124                self.p_func = lambdify((x, y, xi_internal, eta_internal), symbol, 'numpy')
 125            else:
 126                raise ValueError("mode must be 'auto' or 'symbol'")
 127
 128        else:
 129            raise NotImplementedError("Only 1D and 2D supported")
 130
 131        if mode == 'auto':
 132            print("\nsymbol = ")
 133            pprint(self.symbol, num_columns=NUM_COLS)
 134        
 135    def evaluate(self, X, Y, KX, KY, cache=True):
 136        """
 137        Evaluate the pseudo-differential operator's symbol on a grid of spatial and frequency coordinates.
 138
 139        The method dynamically selects between 1D and 2D evaluation based on the spatial dimension.
 140        If caching is enabled and a cached symbol exists, it returns the cached result to avoid recomputation.
 141
 142        Parameters
 143        ----------
 144        X, Y : ndarray
 145            Spatial grid coordinates. In 1D, Y is ignored.
 146        KX, KY : ndarray
 147            Frequency grid coordinates. In 1D, KY is ignored.
 148        cache : bool, default=True
 149            If True, stores the computed symbol for reuse in subsequent calls to avoid redundant computation.
 150
 151        Returns
 152        -------
 153        ndarray
 154            Evaluated symbol values over the input grid. Shape matches the input spatial/frequency grids.
 155
 156        Raises
 157        ------
 158        NotImplementedError
 159            If the spatial dimension is not 1D or 2D.
 160        """
 161        if cache and self.symbol_cached is not None:
 162            return self.symbol_cached
 163
 164        if self.dim == 1:
 165            symbol = self.p_func(X, KX)
 166        elif self.dim == 2:
 167            symbol = self.p_func(X, Y, KX, KY)
 168
 169        if cache:
 170            self.symbol_cached = symbol
 171
 172        return symbol
 173
 174    def clear_cache(self):
 175        """
 176        Clear cached symbol evaluations.
 177        """        
 178        self.symbol_cached = None
 179
 180    def apply(self, u, x_grid, kx, boundary_condition='periodic', 
 181              y_grid=None, ky=None, dealiasing_mask=None,
 182              freq_window='gaussian', clamp=1e6, space_window=False):
 183        """
 184        Apply the pseudo-differential operator to the input field u.
 185    
 186        This method dispatches the application of the pseudo-differential operator based on:
 187        
 188        - Whether the symbol is spatially dependent (x/y)
 189        - The boundary condition in use (periodic or dirichlet)
 190    
 191        Supported operations:
 192        
 193        - Constant-coefficient symbols: applied via Fourier multiplication.
 194        - Spatially varying symbols: applied via Kohn–Nirenberg quantization.
 195        - Dirichlet boundary conditions: handled with non-periodic convolution-like quantization.
 196    
 197        Dispatch Logic:\n
 198        if not self.is_spatial: u ↦ Op(p)(D) ⋅ u = 𝓕⁻¹[ p(ξ) ⋅ 𝓕(u) ]\n
 199        elif periodic: u ↦ Op(p)(x,D) ⋅ u ≈ ∫ eᶦˣᶿ p(x, ξ) 𝓕(u)(ξ) dξ based of FFT (quicker)\n
 200        elif dirichlet: u ↦ Op(p)(x,D) ⋅ u ≈ u ≈ ∫ eᶦˣᶿ p(x, ξ) 𝓕(u)(ξ) dξ (slower)\n
 201        
 202        Parameters
 203        ----------
 204        u : ndarray
 205            Function to which the operator is applied
 206        x_grid : ndarray
 207            Spatial grid in x direction
 208        kx : ndarray
 209            Frequency grid in x direction
 210        boundary_condition : str
 211            'periodic' or 'dirichlet'
 212        y_grid : ndarray, optional
 213            Spatial grid in y direction (for 2D)
 214        ky : ndarray, optional
 215            Frequency grid in y direction (for 2D)
 216        dealiasing_mask : ndarray, optional
 217            Dealiasing mask
 218        freq_window : str
 219            Frequency windowing ('gaussian' or 'hann')
 220        clamp : float
 221            Clamp symbol values to [-clamp, clamp]
 222        space_window : bool
 223            Apply spatial windowing
 224            
 225        Returns
 226        -------
 227        ndarray
 228            Result of applying the operator
 229        """
 230        # Check if symbol depends on spatial variables
 231        is_spatial = self._is_spatial_dependent()
 232        
 233        # Case 1: Constant symbol with periodic BC (fast path)
 234        if not is_spatial and boundary_condition == 'periodic':
 235            return self._apply_constant_fft(u, x_grid, kx, y_grid, ky, dealiasing_mask)
 236        
 237        # Case 2: Spatial symbol with periodic BC
 238        elif boundary_condition == 'periodic':
 239            symbol_func = self._get_symbol_func()
 240            return kohn_nirenberg_fft(
 241                u_vals=u,
 242                symbol_func=symbol_func,
 243                x_grid=x_grid,
 244                kx=kx,
 245                fft_func=self.fft,
 246                ifft_func=self.ifft,
 247                dim=self.dim,
 248                y_grid=y_grid,
 249                ky=ky,
 250                freq_window=freq_window,
 251                clamp=clamp,
 252                space_window=space_window
 253            )
 254        
 255        # Case 3: Dirichlet BC (non-periodic)
 256        elif boundary_condition == 'dirichlet':
 257            symbol_func = self._get_symbol_func()
 258            
 259            if self.dim == 1:
 260                return kohn_nirenberg_nonperiodic(
 261                    u_vals=u,
 262                    x_grid=x_grid,
 263                    xi_grid=kx,
 264                    symbol_func=symbol_func,
 265                    freq_window=freq_window,
 266                    clamp=clamp,
 267                    space_window=space_window
 268                )
 269            elif self.dim == 2:
 270                return kohn_nirenberg_nonperiodic(
 271                    u_vals=u,
 272                    x_grid=(x_grid, y_grid),
 273                    xi_grid=(kx, ky),
 274                    symbol_func=symbol_func,
 275                    freq_window=freq_window,
 276                    clamp=clamp,
 277                    space_window=space_window
 278                )
 279        
 280        else:
 281            raise ValueError(f"Invalid boundary condition '{boundary_condition}'")
 282    
 283    def _is_spatial_dependent(self):
 284        """
 285        Check if the symbol depends on spatial variables.
 286        
 287        Returns
 288        -------
 289        bool
 290            True if symbol depends on x (or x, y)
 291        """
 292        if self.dim == 1:
 293            return self.symbol.has(self.vars_x[0])
 294        elif self.dim == 2:
 295            x, y = self.vars_x
 296            return self.symbol.has(x) or self.symbol.has(y)
 297        else:
 298            return False
 299    
 300    def _get_symbol_func(self):
 301        """
 302        Get a lambdified version of the symbol.
 303        
 304        Returns
 305        -------
 306        callable
 307            Lambdified symbol function
 308        """
 309        if self.dim == 1:
 310            x = self.vars_x[0]
 311            xi = symbols('xi', real=True)
 312            return lambdify((x, xi), self.symbol, 'numpy')
 313        elif self.dim == 2:
 314            x, y = self.vars_x
 315            xi, eta = symbols('xi eta', real=True)
 316            return lambdify((x, y, xi, eta), self.symbol, 'numpy')
 317        else:
 318            raise NotImplementedError("Only 1D and 2D supported")
 319    
 320    def _apply_constant_fft(self, u, x_grid, kx, y_grid, ky, dealiasing_mask):
 321        """
 322        Apply a constant-coefficient pseudo-differential operator in Fourier space.
 323
 324        This method assumes the symbol is diagonal in the Fourier basis and acts as a 
 325        multiplication operator. It performs the operation:
 326        
 327            (ψu)(x) = 𝓕⁻¹[ -σ(k) · 𝓕[u](k) ]
 328
 329        where:
 330        - σ(k) is the combined pseudo-differential operator symbol
 331        - 𝓕 denotes the forward Fourier transform
 332        - 𝓕⁻¹ denotes the inverse Fourier transform
 333
 334        The dealiasing mask is applied before returning to physical space.
 335        
 336        Parameters
 337        ----------
 338        u : ndarray
 339            Input function
 340        x_grid : ndarray
 341            Spatial grid (x)
 342        kx : ndarray
 343            Frequency grid (x)
 344        y_grid : ndarray, optional
 345            Spatial grid (y, for 2D)
 346        ky : ndarray, optional
 347            Frequency grid (y, for 2D)
 348        dealiasing_mask : ndarray, optional
 349            Dealiasing mask
 350            
 351        Returns
 352        -------
 353        ndarray
 354            Result
 355        """
 356        u_hat = self.fft(u)
 357        
 358        # Evaluate symbol at grid points
 359        if self.dim == 1:
 360            X_dummy = np.zeros_like(kx)
 361            symbol_vals = self.p_func(X_dummy, kx)
 362        elif self.dim == 2:
 363            KX, KY = np.meshgrid(kx, ky, indexing='ij')
 364            X_dummy = np.zeros_like(KX)
 365            Y_dummy = np.zeros_like(KY)
 366            symbol_vals = self.p_func(X_dummy, Y_dummy, KX, KY)
 367        else:
 368            raise ValueError("Only 1D and 2D supported")
 369        
 370        # Apply symbol
 371        u_hat *= symbol_vals
 372        
 373        # Apply dealiasing
 374        if dealiasing_mask is not None:
 375            u_hat *= dealiasing_mask
 376        
 377        return self.ifft(u_hat)
 378
 379    def principal_symbol(self, order=1):
 380        """
 381        Compute the leading homogeneous component of the pseudo-differential symbol.
 382
 383        This method extracts the principal part of the symbol, which is the dominant 
 384        term under high-frequency asymptotics (|ξ| → ∞). The expansion is performed 
 385        in polar coordinates for 2D symbols to maintain rotational symmetry, then 
 386        converted back to Cartesian form.
 387
 388        Parameters
 389        ----------
 390        order : int
 391            Order of the asymptotic expansion in powers of 1/ρ, where ρ = |ξ| in 1D 
 392            or ρ = sqrt(ξ² + η²) in 2D. Only the leading-order term is returned.
 393
 394        Returns
 395        -------
 396        sympy.Expr
 397            The principal symbol component, homogeneous of degree `m - order`, where 
 398            `m` is the original symbol's order.
 399
 400        Notes:
 401        - In 1D, uses direct series expansion in ξ.
 402        - In 2D, expands in radial variable ρ while preserving angular dependence.
 403        - Useful for microlocal analysis and constructing parametrices.
 404        """
 405
 406        p = self.symbol
 407        if self.dim == 1:
 408            xi = symbols('xi', real=True, positive=True)
 409            return simplify(series(p, xi, oo, n=order).removeO())
 410        elif self.dim == 2:
 411            xi, eta = symbols('xi eta', real=True, positive=True)
 412            # Homogeneous radial expansion: we set (ξ, η) = ρ (cosθ, sinθ)
 413            rho, theta = symbols('rho theta', real=True, positive=True)
 414            p_rho = p.subs({xi: rho * cos(theta), eta: rho * sin(theta)})
 415            expansion = series(p_rho, rho, oo, n=order).removeO()
 416            # Revert back to (ξ, η)
 417            expansion_cart = expansion.subs({rho: sqrt(xi**2 + eta**2),
 418                                             cos(theta): xi / sqrt(xi**2 + eta**2),
 419                                             sin(theta): eta / sqrt(xi**2 + eta**2)})
 420            return simplify(powdenest(expansion_cart, force=True))
 421                       
 422    def is_homogeneous(self, tol=1e-10):
 423        """
 424        Check whether the symbol is homogeneous in the frequency variables.
 425    
 426        Returns
 427        -------
 428        (bool, Rational or float or None)
 429            Tuple (is_homogeneous, degree) where:
 430            - is_homogeneous: True if the symbol satisfies p(λξ, λη) = λ^m * p(ξ, η)
 431            - degree: the detected degree m if homogeneous, or None
 432        """
 433        from sympy import symbols, simplify, expand, Eq
 434        from sympy.abc import l
 435    
 436        if self.dim == 1:
 437            xi = symbols('xi', real=True, positive=True)
 438            l = symbols('l', real=True, positive=True)
 439            p = self.symbol
 440            p_scaled = p.subs(xi, l * xi)
 441            ratio = simplify(p_scaled / p)
 442            if ratio.has(xi):
 443                return False, None
 444            try:
 445                deg = simplify(ratio).as_base_exp()[1]
 446                return True, deg
 447            except Exception:
 448                return False, None
 449    
 450        elif self.dim == 2:
 451            xi, eta = symbols('xi eta', real=True, positive=True)
 452            l = symbols('l', real=True, positive=True)
 453            p = self.symbol
 454            p_scaled = p.subs({xi: l * xi, eta: l * eta})
 455            ratio = simplify(p_scaled / p)
 456            # If ratio == l**m with no (xi, eta) left, it's homogeneous
 457            if ratio.has(xi, eta):
 458                return False, None
 459            try:
 460                base, exp = ratio.as_base_exp()
 461                if base == l:
 462                    return True, exp
 463            except Exception:
 464                pass
 465            return False, None
 466
 467    def symbol_order(self, max_order=10, tol=1e-3):
 468        """
 469        Estimate the homogeneity order of the pseudo-differential symbol in high-frequency asymptotics.
 470    
 471        This method attempts to determine the leading-order behavior of the symbol p(x, ξ) or p(x, y, ξ, η)
 472        as |ξ| → ∞ (in 1D) or |(ξ, η)| → ∞ (in 2D). The returned value represents the asymptotic growth or decay rate,
 473        which is essential for understanding the regularity and mapping properties of the corresponding operator.
 474    
 475        The function uses symbolic preprocessing to ensure proper factorization of frequency variables,
 476        especially in sqrt and power expressions, to avoid erroneous order detection (e.g., due to hidden scaling).
 477    
 478        Parameters
 479        ----------
 480        max_order : int, optional
 481            Maximum number of terms to consider in the series expansion. Default is 10.
 482        tol : float, optional
 483            Tolerance threshold for evaluating the coefficient magnitude. If the coefficient is too small,
 484            the detected order may be discarded. Default is 1e-3.
 485    
 486        Returns
 487        -------
 488        float or None
 489            - If the symbol is homogeneous, returns its exact homogeneity degree as a float.
 490            - Otherwise, estimates the dominant asymptotic order from leading terms in the expansion.
 491            - Returns None if no valid order could be determined.
 492    
 493        Notes
 494        -----
 495        - In 1D:
 496            Two strategies are used:
 497                1. Expand directly in xi at infinity.
 498                2. Substitute xi = 1/z and expand around z = 0.
 499    
 500        - In 2D:
 501            - Transform the symbol into polar coordinates: (xi, eta) = rho*(cos(theta), sin(theta)).
 502            - Expand in rho at infinity, then extract the leading term's power.
 503            - An alternative substitution using 1/z is also tried if the first method fails.
 504    
 505        - Preprocessing steps:
 506            - Sqrt expressions involving frequencies are rewritten to isolate the leading variable.
 507            - Power expressions are factored explicitly to ensure correct symbolic scaling.
 508    
 509        - If the symbol is not homogeneous, a warning is issued, and the result should be interpreted with care.
 510        
 511        - For non-homogeneous symbols, only the principal asymptotic term is considered.
 512    
 513        Raises
 514        ------
 515        NotImplementedError
 516            If the spatial dimension is neither 1 nor 2.
 517        """
 518        from sympy import (
 519            symbols, series, simplify, sqrt, cos, sin, oo, powdenest, radsimp,
 520            expand, expand_power_base
 521        )
 522    
 523        def preprocess_sqrt(expr, freq):
 524            return expr.replace(
 525                lambda e: e.func == sqrt and freq in e.free_symbols,
 526                lambda e: freq * sqrt(1 + (e.args[0] - freq**2) / freq**2)
 527            )
 528    
 529        def preprocess_power(expr, freq):
 530            return expr.replace(
 531                lambda e: e.is_Pow and freq in e.free_symbols,
 532                lambda e: freq**e.exp * (1 + e.base / freq**e.base.as_powers_dict().get(freq, 0))**e.exp
 533            )
 534    
 535        def validate_order(power, coeff, vars_x, tol):
 536            if power is None:
 537                return None
 538            if any(v in coeff.free_symbols for v in vars_x):
 539                print("⚠️ Coefficient depends on spatial variables; ignoring")
 540                return None
 541            try:
 542                coeff_val = abs(float(coeff.evalf()))
 543                if coeff_val < tol:
 544                    print(f"⚠️ Coefficient too small ({coeff_val:.2e} < {tol})")
 545                    return None
 546            except Exception as e:
 547                print(f"⚠️ Coefficient evaluation failed: {e}")
 548                return None
 549            return int(power) if power == int(power) else float(power)
 550    
 551        # Homogeneity check
 552        is_homog, degree = self.is_homogeneous()
 553        if is_homog:
 554            return float(degree)
 555        else:
 556            print("⚠️ The symbol is not homogeneous. The asymptotic order is not well defined.")
 557    
 558        if self.dim == 1:
 559            x = self.vars_x[0]
 560            xi = symbols('xi', real=True, positive=True)
 561    
 562            try:
 563                print("1D symbol_order - method 1")
 564                expr = preprocess_sqrt(self.symbol, xi)
 565                s = series(expr, xi, oo, n=max_order).removeO()
 566                lead = simplify(powdenest(s.as_leading_term(xi), force=True))
 567                power = lead.as_powers_dict().get(xi, None)
 568                coeff = lead / xi**power if power is not None else 0
 569                print("lead =", lead)
 570                print("power =", power)
 571                print("coeff =", coeff)
 572                order = validate_order(power, coeff, [x], tol)
 573                if order is not None:
 574                    return order
 575            except Exception:
 576                pass
 577    
 578            try:
 579                print("1D symbol_order - method 2")
 580                z = symbols('z', real=True, positive=True)
 581                expr_z = preprocess_sqrt(self.symbol.subs(xi, 1/z), 1/z)
 582                s = series(expr_z, z, 0, n=max_order).removeO()
 583                lead = simplify(powdenest(s.as_leading_term(z), force=True))
 584                power = lead.as_powers_dict().get(z, None)
 585                coeff = lead / z**power if power is not None else 0
 586                print("lead =", lead)
 587                print("power =", power)
 588                print("coeff =", coeff)
 589                order = validate_order(power, coeff, [x], tol)
 590                if order is not None:
 591                    return -order
 592            except Exception as e:
 593                print(f"⚠️ fallback z failed: {e}")
 594            return None
 595    
 596        elif self.dim == 2:
 597            x, y = self.vars_x
 598            xi, eta = symbols('xi eta', real=True, positive=True)
 599            rho, theta = symbols('rho theta', real=True, positive=True)
 600    
 601            try:
 602                print("2D symbol_order - method 1")
 603                p_rho = self.symbol.subs({xi: rho * cos(theta), eta: rho * sin(theta)})
 604                p_rho = preprocess_power(preprocess_sqrt(p_rho, rho), rho)
 605                s = series(simplify(p_rho), rho, oo, n=max_order).removeO()
 606                lead = radsimp(simplify(powdenest(s.as_leading_term(rho), force=True)))
 607                power = lead.as_powers_dict().get(rho, None)
 608                coeff = lead / rho**power if power is not None else 0
 609                print("lead =", lead)
 610                print("power =", power)
 611                print("coeff =", coeff)
 612                order = validate_order(power, coeff, [x, y], tol)
 613                if order is not None:
 614                    return order
 615            except Exception as e:
 616                print(f"⚠️ polar expansion failed: {e}")
 617    
 618            try:
 619                print("2D symbol_order - method 2")
 620                z = symbols('z', real=True, positive=True)
 621                xi_eta = {xi: (1/z) * cos(theta), eta: (1/z) * sin(theta)}
 622                p_rho = preprocess_sqrt(self.symbol.subs(xi_eta), 1/z)
 623                s = series(simplify(p_rho), z, 0, n=max_order).removeO()
 624                lead = radsimp(simplify(powdenest(s.as_leading_term(z), force=True)))
 625                power = lead.as_powers_dict().get(z, None)
 626                coeff = lead / z**power if power is not None else 0
 627                print("lead =", lead)
 628                print("power =", power)
 629                print("coeff =", coeff)
 630                order = validate_order(power, coeff, [x, y], tol)
 631                if order is not None:
 632                    return -order
 633            except Exception as e:
 634                print(f"⚠️ fallback z (2D) failed: {e}")
 635            return None
 636    
 637        else:
 638            raise NotImplementedError("Only 1D and 2D supported.")
 639
 640    
 641    def asymptotic_expansion(self, order=3):
 642        """
 643        Compute the asymptotic expansion of the symbol as |ξ| → ∞ (high-frequency regime).
 644    
 645        This method expands the pseudo-differential symbol in inverse powers of the 
 646        frequency variable(s), either in 1D or 2D. It handles both polynomial and 
 647        exponential symbols by performing a series expansion in 1/|ξ| up to the specified order.
 648    
 649        The expansion is performed directly in Cartesian coordinates for 1D symbols.
 650        For 2D symbols, the method uses polar coordinates (ρ, θ) to perform the expansion 
 651        at infinity in ρ, then converts the result back to Cartesian coordinates.
 652    
 653        Parameters
 654        ----------
 655        order : int, optional
 656            Maximum order of the asymptotic expansion. Default is 3.
 657    
 658        Returns
 659        -------
 660        sympy.Expr
 661            The asymptotic expansion of the symbol up to the given order, expressed in Cartesian coordinates.
 662            If expansion fails, returns the original unexpanded symbol.
 663    
 664        Notes:
 665        - In 1D: expansion is performed directly in terms of ξ.
 666        - In 2D: the symbol is first rewritten in polar coordinates (ρ,θ), expanded asymptotically 
 667          in ρ → ∞, then converted back to Cartesian coordinates (ξ,η).
 668        - Handles special case when the symbol is an exponential function by expanding its argument.
 669        - Symbolic normalization is applied early (via `simplify`) for 2D expressions to improve convergence.
 670        - Robust to failures: catches exceptions and issues warnings instead of raising errors.
 671        - Final expression is simplified using `powdenest` and `expand` for improved readability.
 672        """
 673        p = self.symbol
 674    
 675        if self.dim == 1:
 676            xi = symbols('xi', real=True, positive=True)
 677    
 678            try:
 679                # Case: exponential function
 680                if p.func == exp and len(p.args) == 1:
 681                    arg = p.args[0]
 682                    arg_series = series(arg, xi, oo, n=order).removeO()
 683                    expanded = series(exp(expand(arg_series)), xi, oo, n=order).removeO()
 684                    return simplify(powdenest(expanded, force=True))
 685                else:
 686                    expanded = series(p, xi, oo, n=order).removeO()
 687                    return simplify(powdenest(expanded, force=True))
 688    
 689            except Exception as e:
 690                print(f"Warning: 1D expansion failed: {e}")
 691                return p
 692    
 693        elif self.dim == 2:
 694            xi, eta = symbols('xi eta', real=True, positive=True)
 695            rho, theta = symbols('rho theta', real=True, positive=True)
 696    
 697            # Normalize before substitution
 698            p = simplify(p)
 699    
 700            # Substitute polar coordinates
 701            p_polar = p.subs({
 702                xi: rho * cos(theta),
 703                eta: rho * sin(theta)
 704            })
 705    
 706            try:
 707                # Handle exponentials
 708                if p_polar.func == exp and len(p_polar.args) == 1:
 709                    arg = p_polar.args[0]
 710                    arg_series = series(arg, rho, oo, n=order).removeO()
 711                    expanded = series(exp(expand(arg_series)), rho, oo, n=order).removeO()
 712                else:
 713                    expanded = series(p_polar, rho, oo, n=order).removeO()
 714    
 715                # Convert back to Cartesian
 716                norm = sqrt(xi**2 + eta**2)
 717                expansion_cart = expanded.subs({
 718                    rho: norm,
 719                    cos(theta): xi / norm,
 720                    sin(theta): eta / norm
 721                })
 722    
 723                # Final simplifications
 724                result = simplify(powdenest(expansion_cart, force=True))
 725                result = expand(result)
 726                return result
 727    
 728            except Exception as e:
 729                print(f"Warning: 2D expansion failed: {e}")
 730                return p  
 731            
 732    def compose_asymptotic(self, other, order=1, mode='kn', sign_convention=None):
 733        """
 734        Compose two pseudo-differential operators using an asymptotic expansion
 735        in the chosen quantization scheme (Kohn–Nirenberg or Weyl).
 736    
 737        Parameters
 738        ----------
 739        other : PseudoDifferentialOperator
 740            The operator to compose with this one.
 741        order : int, default=1
 742            Maximum order of the asymptotic expansion.
 743        mode : {'kn', 'weyl'}, default='kn'
 744            Quantization mode:
 745            - 'kn' : Kohn–Nirenberg quantization (left-quantized)
 746            - 'weyl' : Weyl symmetric quantization
 747        sign_convention : {'standard', 'inverse'}, optional
 748            Controls the phase factor convention for the KN case:
 749            - 'standard' → (i)^(-n), gives [x, ξ] = +i (physics convention)
 750            - 'inverse' → (i)^(+n), gives [x, ξ] = -i (mathematical adjoint convention)
 751            If None, defaults to 'standard'.
 752    
 753        Returns
 754        -------
 755        sympy.Expr
 756            Symbolic expression for the composed symbol up to the given order.
 757    
 758        Notes
 759        -----
 760        - In 1D (Kohn–Nirenberg):
 761            (p ∘ q)(x, ξ) ~ Σₙ (1/n!) (i sgn)^n ∂_ξⁿ p(x, ξ) ∂_xⁿ q(x, ξ)
 762        - In 1D (Weyl):
 763            (p # q)(x, ξ) = exp[(i/2)(∂_ξ^p ∂_x^q - ∂_x^p ∂_ξ^q)] p(x, ξ) q(x, ξ)
 764            truncated at given order.
 765    
 766        Examples
 767        --------
 768        X = a*x, Y = b*ξ
 769        X_op.compose_asymptotic(Y_op, order=3, mode='weyl')
 770        """
 771    
 772        from sympy import diff, factorial, simplify, symbols
 773    
 774        assert self.dim == other.dim, "Operator dimensions must match"
 775        p, q = self.symbol, other.symbol
 776    
 777        # Default sign convention
 778        if sign_convention is None:
 779            sign_convention = 'standard'
 780        sign = -1 if sign_convention == 'standard' else +1
 781    
 782        # --- 1D case ---
 783        if self.dim == 1:
 784            x = self.vars_x[0]
 785            xi = symbols('xi', real=True)
 786            result = 0
 787    
 788            if mode == 'kn':  # Kohn–Nirenberg
 789                for n in range(order + 1):
 790                    term = (1 / factorial(n)) * diff(p, xi, n) * diff(q, x, n) * (1j) ** (sign * n)
 791                    result += term
 792    
 793            elif mode == 'weyl':  # Weyl symmetric composition
 794                # Weyl star product: exp((i/2)(∂_ξ^p ∂_x^q - ∂_x^p ∂_ξ^q))
 795                result = 0
 796                for n in range(order + 1):
 797                    for k in range(n + 1):
 798                        # k derivatives acting as (∂_ξ^k p)(∂_x^(n−k) q)
 799                        coeff = (1 / (factorial(k) * factorial(n - k))) * ((1j / 2) ** n) * ((-1) ** (n - k))
 800                        term = coeff * diff(p, xi, k, x, n - k, evaluate=True) * diff(q, x, k, xi, n - k, evaluate=True)
 801                        result += term
 802    
 803            else:
 804                raise ValueError("mode must be either 'kn' or 'weyl'")
 805    
 806            return simplify(result)
 807    
 808        # --- 2D case ---
 809        elif self.dim == 2:
 810            x, y = self.vars_x
 811            xi, eta = symbols('xi eta', real=True)
 812            result = 0
 813    
 814            if mode == 'kn':
 815                for n in range(order + 1):
 816                    for i in range(n + 1):
 817                        j = n - i
 818                        term = (1 / (factorial(i) * factorial(j))) * \
 819                               diff(p, xi, i, eta, j) * diff(q, x, i, y, j) * (1j) ** (sign * n)
 820                        result += term
 821    
 822            elif mode == 'weyl':
 823                for n in range(order + 1):
 824                    for i in range(n + 1):
 825                        j = n - i
 826                        coeff = (1 / (factorial(i) * factorial(j))) * ((1j / 2) ** n) * ((-1) ** (n - i))
 827                        term = coeff * diff(p, xi, i, eta, j, x, 0, y, 0) * diff(q, x, i, y, j, xi, 0, eta, 0)
 828                        result += term
 829            else:
 830                raise ValueError("mode must be either 'kn' or 'weyl'")
 831    
 832            return simplify(result)
 833    
 834        else:
 835            raise NotImplementedError("Only 1D and 2D cases are implemented")
 836
 837    def commutator_symbolic(self, other, order=1, mode='kn', sign_convention=None):
 838        """
 839        Compute the symbolic commutator [A, B] = A∘B − B∘A of two pseudo-differential operators
 840        using formal asymptotic expansion of their composition symbols.
 841    
 842        This method computes the asymptotic expansion of the commutator's symbol up to a given 
 843        order, based on the symbolic calculus of pseudo-differential operators in the 
 844        Kohn–Nirenberg quantization. The result is a purely symbolic sympy expression that 
 845        captures the leading-order noncommutativity of the operators.
 846    
 847        Parameters
 848        ----------
 849        other : PseudoDifferentialOperator
 850            The pseudo-differential operator B to commute with this operator A.
 851        order : int, default=1
 852            Maximum order of the asymptotic expansion. 
 853            - order=1 yields the leading term proportional to the Poisson bracket {p, q}.
 854            - Higher orders include correction terms involving higher mixed derivatives.
 855    
 856        Returns
 857        -------
 858        sympy.Expr
 859            Symbolic expression for the asymptotic expansion of the commutator symbol 
 860            σ([A,B]) = σ(A∘B − B∘A).
 861    
 862        """
 863        assert self.dim == other.dim, "Operator dimensions must match"
 864        p, q = self.symbol, other.symbol
 865    
 866        pq = self.compose_asymptotic(other, order=order, mode=mode, sign_convention=sign_convention)
 867        qp = other.compose_asymptotic(self, order=order, mode=mode, sign_convention=sign_convention)
 868        
 869        comm_symbol = simplify(pq-qp)
 870
 871        return comm_symbol
 872
 873    def right_inverse_asymptotic(self, order=1):
 874        """
 875        Construct a formal right inverse R of the pseudo-differential operator P such that 
 876        the composition P ∘ R equals the identity plus a smoothing operator of order -order.
 877    
 878        This method computes an asymptotic expansion for the right inverse using recursive 
 879        corrections based on derivatives of the symbol p(x, ξ) and lower-order terms of R.
 880    
 881        Parameters
 882        ----------
 883        order : int
 884            Number of terms to include in the asymptotic expansion. Higher values improve 
 885            approximation at the cost of complexity and computational effort.
 886    
 887        Returns
 888        -------
 889        sympy.Expr
 890            The symbolic expression representing the formal right inverse R(x, ξ), which satisfies:
 891            P ∘ R = Id + O(⟨ξ⟩^{-order}), where ⟨ξ⟩ = (1 + |ξ|²)^{1/2}.
 892    
 893        Notes
 894        -----
 895        - In 1D: The recursion involves spatial derivatives of R and derivatives of p with respect to ξ.
 896        - In 2D: The multi-index generalization is used with mixed derivatives in ξ and η.
 897        - The construction relies on the non-vanishing of the principal symbol p to ensure invertibility.
 898        - Each term in the expansion corresponds to higher-order corrections involving commutators 
 899          between the operator P and the current approximation of R.
 900        """
 901        p = self.symbol
 902        if self.dim == 1:
 903            x = self.vars_x[0]
 904            xi = symbols('xi', real=True)
 905            r = 1 / p.subs(xi, xi)  # r0
 906            R = r
 907            for n in range(1, order + 1):
 908                term = 0
 909                for k in range(1, n + 1):
 910                    coeff = (1j)**(-k) / factorial(k)
 911                    inner = diff(p, xi, k) * diff(R, x, k)
 912                    term += coeff * inner
 913                R = R - r * term
 914        elif self.dim == 2:
 915            x, y = self.vars_x
 916            xi, eta = symbols('xi eta', real=True)
 917            r = 1 / p.subs({xi: xi, eta: eta})
 918            R = r
 919            for n in range(1, order + 1):
 920                term = 0
 921                for k1 in range(n + 1):
 922                    for k2 in range(n + 1 - k1):
 923                        if k1 + k2 == 0: continue
 924                        coeff = (1j)**(-(k1 + k2)) / (factorial(k1) * factorial(k2))
 925                        dp = diff(p, xi, k1, eta, k2)
 926                        dR = diff(R, x, k1, y, k2)
 927                        term += coeff * dp * dR
 928                R = R - r * term
 929        return R
 930
 931    def left_inverse_asymptotic(self, order=1):
 932        """
 933        Construct a formal left inverse L such that the composition L ∘ P equals the identity 
 934        operator up to terms of order ξ^{-order}. This expansion is performed asymptotically 
 935        at infinity in the frequency variable(s).
 936    
 937        The left inverse is built iteratively using symbolic differentiation and the 
 938        method of asymptotic expansions for pseudo-differential operators. It ensures that:
 939        
 940            L(P(x,ξ),x,D) ∘ P(x,D) = Id + smoothing operator of order -order
 941    
 942        Parameters
 943        ----------
 944        order : int, optional
 945            Maximum number of terms in the asymptotic expansion (default is 1). Higher values 
 946            yield more accurate inverses at the cost of increased computational complexity.
 947    
 948        Returns
 949        -------
 950        sympy.Expr
 951            Symbolic expression representing the principal symbol of the formal left inverse 
 952            operator L(x,ξ). This expression depends on spatial variables and frequencies, 
 953            and includes correction terms up to the specified order.
 954    
 955        Notes
 956        -----
 957        - In 1D: Uses recursive application of the Leibniz formula for symbols.
 958        - In 2D: Generalizes to multi-indices for mixed derivatives in (x,y) and (ξ,η).
 959        - Each term involves combinations of derivatives of the original symbol p(x,ξ) and 
 960          previously computed terms of the inverse.
 961        - Coefficients include powers of 1j (i) and factorial normalization for derivative terms.
 962        """
 963        p = self.symbol
 964        if self.dim == 1:
 965            x = self.vars_x[0]
 966            xi = symbols('xi', real=True)
 967            l = 1 / p.subs(xi, xi)
 968            L = l
 969            for n in range(1, order + 1):
 970                term = 0
 971                for k in range(1, n + 1):
 972                    coeff = (1j)**(-k) / factorial(k)
 973                    inner = diff(L, xi, k) * diff(p, x, k)
 974                    term += coeff * inner
 975                L = L - term * l
 976        elif self.dim == 2:
 977            x, y = self.vars_x
 978            xi, eta = symbols('xi eta', real=True)
 979            l = 1 / p.subs({xi: xi, eta: eta})
 980            L = l
 981            for n in range(1, order + 1):
 982                term = 0
 983                for k1 in range(n + 1):
 984                    for k2 in range(n + 1 - k1):
 985                        if k1 + k2 == 0: continue
 986                        coeff = (1j)**(-(k1 + k2)) / (factorial(k1) * factorial(k2))
 987                        dp = diff(p, x, k1, y, k2)
 988                        dL = diff(L, xi, k1, eta, k2)
 989                        term += coeff * dL * dp
 990                L = L - term * l
 991        return L
 992
 993    def formal_adjoint(self):
 994        """
 995        Compute the formal adjoint symbol P* of the pseudo-differential operator.
 996
 997        The adjoint is defined such that for any test functions u and v,
 998        ⟨P u, v⟩ = ⟨u, P* v⟩ holds in the distributional sense. This is obtained by 
 999        taking the complex conjugate of the symbol and expanding it asymptotically 
1000        at infinity to ensure proper behavior under integration by parts.
1001
1002        Returns
1003        -------
1004        sympy.Expr
1005            The adjoint symbol P*(x, ξ) in 1D or P*(x, y, ξ, η) in 2D.
1006        
1007        Notes:
1008        - In 1D, the expansion is performed in powers of 1/|ξ|.
1009        - In 2D, the expansion is radial in |ξ| = sqrt(ξ² + η²).
1010        - This method ensures symbolic simplifications for readability and efficiency.
1011        """
1012        p = self.symbol
1013        if self.dim == 1:
1014            x, = self.vars_x
1015            xi = symbols('xi', real=True)
1016            p_star = conjugate(p)
1017            p_star = simplify(series(p_star, xi, oo, n=6).removeO())
1018            return p_star
1019        elif self.dim == 2:
1020            x, y = self.vars_x
1021            xi, eta = symbols('xi eta', real=True)
1022            p_star = conjugate(p)
1023            p_star = simplify(series(p_star, sqrt(xi**2 + eta**2), oo, n=6).removeO())
1024            return p_star
1025
1026    def exponential_symbol(self, t=1.0, order=1, mode='kn', sign_convention=None):
1027        """
1028        Compute the symbol of exp(tP) using asymptotic expansion methods.
1029        
1030        This method calculates the exponential of a pseudo-differential operator 
1031        using either a direct power series expansion or a Magnus expansion, 
1032        depending on the structure of the symbol. The result is valid up to 
1033        the specified asymptotic order.
1034        
1035        Parameters
1036        ----------
1037        t : float or sympy.Symbol, default=1.0
1038            Time or evolution parameter. Common uses:
1039            - t = -i*τ for Schrödinger evolution: exp(-iτH)
1040            - t = τ for heat/diffusion: exp(τΔ)
1041            - t for general propagators
1042        order : int, default=3
1043            Maximum order of the asymptotic expansion. Higher orders include 
1044            more composition terms, improving accuracy for small t or when 
1045            non-commutativity effects are significant.
1046        
1047        Returns
1048        -------
1049        sympy.Expr
1050            Symbolic expression for the exponential operator symbol, computed 
1051            as an asymptotic series up to the specified order.
1052        
1053        Notes
1054        -----
1055        - For commutative symbols (e.g., pure multiplication operators), the 
1056          exponential is exact: exp(tP) = exp(t*p(x,ξ)).
1057        
1058        - For general non-commutative operators, the method uses the BCH-type 
1059          expansion via iterated composition:
1060          exp(tP) ~ I + tP + (t²/2!)P∘P + (t³/3!)P∘P∘P + ...
1061          
1062        - Each power P^n is computed via compose_asymptotic, which accounts 
1063          for the non-commutativity through derivative terms.
1064        
1065        - The expansion is valid for |t| small enough or when the symbol has 
1066          appropriate decay/growth properties.
1067        
1068        - In quantum mechanics (Schrödinger): U(t) = exp(-itH/ℏ) represents 
1069          the time evolution operator.
1070        
1071        - In parabolic PDEs (heat equation): exp(tΔ) is the heat kernel.
1072
1073        """
1074        if self.dim == 1:
1075            x = self.vars_x[0]
1076            xi = symbols('xi', real=True)
1077            
1078            # Initialize with identity
1079            result = 1
1080            
1081            # First order term: tP
1082            current_power = self.symbol
1083            result += t * current_power
1084            
1085            # Higher order terms: (t^n/n!) P^n computed via composition
1086            for n in range(2, order + 1):
1087                # Compute P^n = P^(n-1) ∘ P via asymptotic composition
1088                # We use a temporary operator for composition
1089                temp_op = PseudoDifferentialOperator(
1090                    current_power, [x], mode='symbol'
1091                )
1092                current_power = temp_op.compose_asymptotic(self, order=order, mode=mode, sign_convention=sign_convention)
1093                
1094                # Add term (t^n/n!) * P^n
1095                coeff = t**n / factorial(n)
1096                result += coeff * current_power
1097            
1098            return simplify(result)
1099        
1100        elif self.dim == 2:
1101            x, y = self.vars_x
1102            xi, eta = symbols('xi eta', real=True)
1103            
1104            # Initialize with identity
1105            result = 1
1106            
1107            # First order term: tP
1108            current_power = self.symbol
1109            result += t * current_power
1110            
1111            # Higher order terms: (t^n/n!) P^n computed via composition
1112            for n in range(2, order + 1):
1113                # Compute P^n = P^(n-1) ∘ P via asymptotic composition
1114                temp_op = PseudoDifferentialOperator(
1115                    current_power, [x, y], mode='symbol'
1116                )
1117                current_power = temp_op.compose_asymptotic(self, order=order, mode=mode, sign_convention=sign_convention)
1118                
1119                # Add term (t^n/n!) * P^n
1120                coeff = t**n / factorial(n)
1121                result += coeff * current_power
1122            
1123            return simplify(result)
1124        
1125        else:
1126            raise NotImplementedError("Only 1D and 2D operators are supported")
1127        
1128    def trace_formula(self, volume_element=None, numerical=False, 
1129                      x_bounds=None, xi_bounds=None):
1130        """
1131        Compute the semiclassical trace of the pseudo-differential operator.
1132        
1133        The trace formula relates the quantum trace of an operator to a 
1134        phase-space integral of its symbol, providing a fundamental link 
1135        between classical and quantum mechanics. This implementation supports 
1136        both symbolic and numerical integration.
1137        
1138        Parameters
1139        ----------
1140        volume_element : sympy.Expr, optional
1141            Custom volume element for the phase space integration. If None, 
1142            uses the standard Liouville measure dx dξ/(2π)^d.
1143        numerical : bool, default=False
1144            If True, perform numerical integration over specified bounds.
1145            If False, attempt symbolic integration (may fail for complex symbols).
1146        x_bounds : tuple of tuples, optional
1147            Spatial integration bounds. For 1D: ((x_min, x_max),)
1148            For 2D: ((x_min, x_max), (y_min, y_max))
1149            Required if numerical=True.
1150        xi_bounds : tuple of tuples, optional
1151            Frequency integration bounds. For 1D: ((xi_min, xi_max),)
1152            For 2D: ((xi_min, xi_max), (eta_min, eta_max))
1153            Required if numerical=True.
1154        
1155        Returns
1156        -------
1157        sympy.Expr or float
1158            The trace of the operator. Returns a symbolic expression if 
1159            numerical=False, or a float if numerical=True.
1160        
1161        Notes
1162        -----
1163        - The semiclassical trace formula states:
1164          Tr(P) = (2π)^{-d} ∫∫ p(x,ξ) dx dξ
1165          where d is the spatial dimension and p(x,ξ) is the operator symbol.
1166        
1167        - For 1D: Tr(P) = (1/2π) ∫_{-∞}^{∞} ∫_{-∞}^{∞} p(x,ξ) dx dξ
1168        
1169        - For 2D: Tr(P) = (1/4π²) ∫∫∫∫ p(x,y,ξ,η) dx dy dξ dη
1170        
1171        - This formula is exact for trace-class operators and provides an 
1172          asymptotic approximation for general pseudo-differential operators.
1173        
1174        - Physical interpretation: the trace counts the "number of states" 
1175          weighted by the observable p(x,ξ).
1176        
1177        - For projection operators (χ_Ω with χ² = χ), the trace gives the 
1178          dimension of the range, related to the phase space volume of Ω.
1179        
1180        - The factor (2π)^{-d} comes from the quantum normalization of 
1181          coherent states / Weyl quantization.
1182        """
1183        from sympy import integrate, simplify, lambdify
1184        from scipy.integrate import dblquad, nquad
1185        
1186        p = self.symbol
1187        
1188        if numerical:
1189            if x_bounds is None or xi_bounds is None:
1190                raise ValueError(
1191                    "x_bounds and xi_bounds must be provided for numerical integration"
1192                )
1193        
1194        if self.dim == 1:
1195            x, = self.vars_x
1196            xi = symbols('xi', real=True)
1197            
1198            if volume_element is None:
1199                volume_element = 1 / (2 * pi)
1200            
1201            if numerical:
1202                # Numerical integration
1203                p_func = lambdify((x, xi), p, 'numpy')
1204                (x_min, x_max), = x_bounds
1205                (xi_min, xi_max), = xi_bounds
1206                
1207                def integrand(xi_val, x_val):
1208                    return p_func(x_val, xi_val)
1209                
1210                result, error = dblquad(
1211                    integrand,
1212                    x_min, x_max,
1213                    lambda x: xi_min, lambda x: xi_max
1214                )
1215                
1216                result *= float(volume_element)
1217                print(f"Numerical trace = {result:.6e} ± {error:.6e}")
1218                return result
1219            
1220            else:
1221                # Symbolic integration
1222                integrand = p * volume_element
1223                
1224                try:
1225                    # Try to integrate over xi first, then x
1226                    integral_xi = integrate(integrand, (xi, -oo, oo))
1227                    integral_x = integrate(integral_xi, (x, -oo, oo))
1228                    return simplify(integral_x)
1229                except:
1230                    print("Warning: Symbolic integration failed. Try numerical=True")
1231                    return integrate(integrand, (xi, -oo, oo), (x, -oo, oo))
1232        
1233        elif self.dim == 2:
1234            x, y = self.vars_x
1235            xi, eta = symbols('xi eta', real=True)
1236            
1237            if volume_element is None:
1238                volume_element = 1 / (4 * pi**2)
1239            
1240            if numerical:
1241                # Numerical integration in 4D
1242                p_func = lambdify((x, y, xi, eta), p, 'numpy')
1243                (x_min, x_max), (y_min, y_max) = x_bounds
1244                (xi_min, xi_max), (eta_min, eta_max) = xi_bounds
1245                
1246                def integrand(eta_val, xi_val, y_val, x_val):
1247                    return p_func(x_val, y_val, xi_val, eta_val)
1248                
1249                result, error = nquad(
1250                    integrand,
1251                    [
1252                        [eta_min, eta_max],
1253                        [xi_min, xi_max],
1254                        [y_min, y_max],
1255                        [x_min, x_max]
1256                    ]
1257                )
1258                
1259                result *= float(volume_element)
1260                print(f"Numerical trace = {result:.6e} ± {error:.6e}")
1261                return result
1262            
1263            else:
1264                # Symbolic integration
1265                integrand = p * volume_element
1266                
1267                try:
1268                    # Integrate in order: eta, xi, y, x
1269                    integral_eta = integrate(integrand, (eta, -oo, oo))
1270                    integral_xi = integrate(integral_eta, (xi, -oo, oo))
1271                    integral_y = integrate(integral_xi, (y, -oo, oo))
1272                    integral_x = integrate(integral_y, (x, -oo, oo))
1273                    return simplify(integral_x)
1274                except:
1275                    print("Warning: Symbolic integration failed. Try numerical=True")
1276                    return integrate(
1277                        integrand,
1278                        (eta, -oo, oo), (xi, -oo, oo),
1279                        (y, -oo, oo), (x, -oo, oo)
1280                    )
1281        
1282        else:
1283            raise NotImplementedError("Only 1D and 2D operators are supported")
1284
1285    def pseudospectrum_analysis(self, x_grid, lambda_real_range, lambda_imag_range,
1286                               epsilon_levels=[0.1, 0.01, 0.001, 0.0001],
1287                               resolution=100, method='spectral', L=None, N=None,
1288                               use_sparse=False, parallel=True, n_workers=4,
1289                               adaptive=False, adaptive_threshold=0.5,
1290                               auto_range=True, plot=True):
1291        """
1292        Compute and visualize the pseudospectrum of the operator.
1293        
1294        Optimizations:
1295        - Uses apply() method instead of manual loops
1296        - Parallel computation of resolvent norms
1297        - Sparse matrix support for large N
1298        - Optional adaptive grid refinement
1299        
1300        Parameters
1301        ----------
1302        x_grid : array
1303            Spatial grid for quantization
1304        lambda_real_range : tuple
1305            (min, max) for real part of λ
1306        lambda_imag_range : tuple
1307            (min, max) for imaginary part of λ
1308        epsilon_levels : list
1309            Levels for ε-pseudospectrum contours
1310        resolution : int
1311            Grid resolution for λ sampling
1312        method : str
1313            'spectral' or 'finite_difference'
1314        L : float, optional
1315            Domain half-length for spectral method
1316        N : int, optional
1317            Number of grid points
1318        use_sparse : bool
1319            Use sparse matrices for large N
1320        parallel : bool
1321            Enable parallel computation
1322        n_workers : int
1323            Number of parallel workers
1324        adaptive : bool
1325            Use adaptive grid refinement
1326        adaptive_threshold : float
1327            Threshold for adaptive refinement
1328            
1329        Returns
1330        -------
1331        dict
1332            Dictionary with pseudospectrum data and operator matrix
1333        """
1334        if self.dim != 1:
1335            raise NotImplementedError('Pseudospectrum analysis currently supports 1D only')
1336        
1337        # Step 1: Build operator matrix
1338        print(f"Building operator matrix using '{method}' method...")
1339        H, x_grid_used, k_grid = self._build_operator_matrix(x_grid, method, L, N)
1340        N_actual = H.shape[0]
1341        
1342        # Step 1.5: Compute eigenvalues FIRST to adjust range if needed
1343        print('Computing eigenvalues...')
1344        eigenvalues = self._compute_eigenvalues(H, use_sparse)
1345        
1346        # Auto-adjust range if requested
1347        if auto_range and eigenvalues is not None:
1348            eig_real_min, eig_real_max = eigenvalues.real.min(), eigenvalues.real.max()
1349            eig_imag_min, eig_imag_max = eigenvalues.imag.min(), eigenvalues.imag.max()
1350            
1351            # Add 20% margin around eigenvalues
1352            margin_real = 0.2 * (eig_real_max - eig_real_min + 1)
1353            margin_imag = max(0.2 * (eig_imag_max - eig_imag_min + 1), 2.0)
1354            
1355            lambda_real_range = (eig_real_min - margin_real, eig_real_max + margin_real)
1356            lambda_imag_range = (eig_imag_min - margin_imag, eig_imag_max + margin_imag)
1357            
1358            print(f'Auto-adjusted λ range:')
1359            print(f'  Re(λ) ∈ [{lambda_real_range[0]:.2f}, {lambda_real_range[1]:.2f}]')
1360            print(f'  Im(λ) ∈ [{lambda_imag_range[0]:.2f}, {lambda_imag_range[1]:.2f}]')
1361        
1362        # Step 2: Compute pseudospectrum with corrected range
1363        print(f'Computing pseudospectrum over {resolution}×{resolution} grid...')
1364        if adaptive:
1365            print('Using adaptive grid refinement...')
1366            Lambda, resolvent_norm, sigma_min_grid = self._compute_pseudospectrum_adaptive(
1367                H, lambda_real_range, lambda_imag_range, resolution,
1368                use_sparse=use_sparse, parallel=parallel, n_workers=n_workers,
1369                threshold=adaptive_threshold
1370            )
1371        else:
1372            Lambda, resolvent_norm, sigma_min_grid = self._compute_pseudospectrum(
1373                H, lambda_real_range, lambda_imag_range, resolution,
1374                use_sparse=use_sparse, parallel=parallel, n_workers=n_workers
1375            )
1376        
1377        # Step 3: Visualize
1378        if plot:
1379            self._plot_pseudospectrum(Lambda, resolvent_norm, sigma_min_grid,
1380                                      epsilon_levels, eigenvalues)
1381        
1382        return {
1383            'lambda_grid': Lambda,
1384            'resolvent_norm': resolvent_norm,
1385            'sigma_min': sigma_min_grid,
1386            'epsilon_levels': epsilon_levels,
1387            'eigenvalues': eigenvalues,
1388            'operator_matrix': H,
1389            'x_grid': x_grid_used,
1390            'k_grid': k_grid
1391        }
1392
1393
1394    def _build_operator_matrix(self, x_grid, method, L, N):
1395        """
1396        Build the discrete operator matrix H.
1397        
1398        Optimized to use the apply() method instead of manual integration.
1399        
1400        Parameters
1401        ----------
1402        x_grid : array
1403            Input spatial grid
1404        method : str
1405            'spectral' or 'finite_difference'
1406        L : float, optional
1407            Domain half-length
1408        N : int, optional
1409            Number of grid points
1410            
1411        Returns
1412        -------
1413        H : ndarray
1414            Operator matrix (N×N)
1415        x_grid_used : ndarray
1416            Actual spatial grid used
1417        k_grid : ndarray
1418            Frequency grid
1419        """
1420        if method == 'spectral':
1421            # Setup spectral grid
1422            if L is None:
1423                L = (x_grid[-1] - x_grid[0]) / 2.0
1424            if N is None:
1425                N = len(x_grid)
1426            
1427            x_grid_spectral = np.linspace(-L, L, N, endpoint=False)
1428            dx = x_grid_spectral[1] - x_grid_spectral[0]
1429            k = np.fft.fftfreq(N, d=dx) * 2.0 * np.pi
1430            
1431            # Build matrix by applying operator to canonical basis
1432            # This is the KEY OPTIMIZATION: use apply() instead of manual loops
1433            H = np.zeros((N, N), dtype=complex)
1434            
1435            for j in range(N):
1436                # Create basis vector e_j
1437                e_j = np.zeros(N, dtype=complex)
1438                e_j[j] = 1.0
1439                
1440                # Apply operator using the existing apply() method
1441                # This automatically handles the symbol evaluation and FFT operations
1442                H[:, j] = self.apply(
1443                    e_j, 
1444                    x_grid_spectral, 
1445                    k,
1446                    boundary_condition='periodic'
1447                )
1448            
1449            print(f'Operator quantized via apply() method: {N}×{N} matrix')
1450            return H, x_grid_spectral, k
1451            
1452        elif method == 'finite_difference':
1453            # Fallback to finite difference (keep original implementation)
1454            N = len(x_grid)
1455            dx = x_grid[1] - x_grid[0]
1456            H = np.zeros((N, N), dtype=complex)
1457            
1458            for i in range(N):
1459                for j in range(N):
1460                    if i == j:
1461                        H[i, j] = self.p_func(x_grid[i], 0.0)
1462                    elif abs(i - j) == 1:
1463                        xi_approx = np.pi / dx
1464                        H[i, j] = self.p_func(
1465                            (x_grid[i] + x_grid[j]) / 2,
1466                            xi_approx * np.sign(i - j)
1467                        ) / (2 * dx)
1468                    elif abs(i - j) == 2:
1469                        xi_approx = 2 * np.pi / dx
1470                        H[i, j] = self.p_func(
1471                            (x_grid[i] + x_grid[j]) / 2,
1472                            xi_approx
1473                        ) / dx ** 2
1474            
1475            print(f'Operator quantized via finite differences: {N}×{N} matrix')
1476            k = np.fft.fftfreq(N, d=dx) * 2.0 * np.pi
1477            return H, x_grid, k
1478            
1479        else:
1480            raise ValueError("method must be 'spectral' or 'finite_difference'")
1481
1482    def _compute_pseudospectrum(self, H, lambda_real_range, lambda_imag_range,
1483                               resolution, use_sparse=False, parallel=True,
1484                               n_workers=4):
1485        """
1486        Compute pseudospectrum on a uniform grid.
1487        
1488        Optimized with parallel computation and optional sparse matrices.
1489        
1490        Parameters
1491        ----------
1492        H : ndarray or sparse matrix
1493            Operator matrix
1494        lambda_real_range : tuple
1495            Range for Re(λ)
1496        lambda_imag_range : tuple
1497            Range for Im(λ)
1498        resolution : int
1499            Grid resolution
1500        use_sparse : bool
1501            Use sparse SVD for large matrices
1502        parallel : bool
1503            Enable parallel computation
1504        n_workers : int
1505            Number of parallel workers
1506            
1507        Returns
1508        -------
1509        Lambda : ndarray
1510            Complex grid of λ values
1511        resolvent_norm : ndarray
1512            Norm of (H - λI)^{-1}
1513        sigma_min_grid : ndarray
1514            Smallest singular value σ_min(H - λI)
1515        """
1516        from scipy.linalg import svdvals
1517        
1518        N = H.shape[0]
1519        lambda_re = np.linspace(*lambda_real_range, resolution)
1520        lambda_im = np.linspace(*lambda_imag_range, resolution)
1521        Lambda_re, Lambda_im = np.meshgrid(lambda_re, lambda_im)
1522        Lambda = Lambda_re + 1j * Lambda_im
1523        
1524        resolvent_norm = np.zeros_like(Lambda, dtype=float)
1525        sigma_min_grid = np.zeros_like(Lambda, dtype=float)
1526        
1527        I = np.eye(N)
1528        
1529        # Convert to sparse if requested and beneficial
1530        if use_sparse and N > 100:
1531            from scipy.sparse import csr_matrix, eye as sparse_eye
1532            from scipy.sparse.linalg import svds
1533            H_sparse = csr_matrix(H)
1534            I_sparse = sparse_eye(N, format='csr')
1535            use_sparse_svd = True
1536            print(f'Using sparse matrices (N={N})')
1537        else:
1538            use_sparse_svd = False
1539        
1540        if parallel and resolution * resolution > 100:
1541            # Parallel computation
1542            Lambda_flat = Lambda.ravel()
1543            
1544            def compute_single_point(idx):
1545                """Compute resolvent norm for a single λ value"""
1546                lam = Lambda_flat[idx]
1547                try:
1548                    if use_sparse_svd:
1549                        # Sparse SVD: compute only smallest singular value
1550                        A = H_sparse - lam * I_sparse
1551                        try:
1552                            # svds can be unstable, wrap in try-except
1553                            s_min = svds(A, k=1, which='SM', 
1554                                       return_singular_vectors=False)[0]
1555                        except:
1556                            # Fallback to dense computation
1557                            s = svdvals(A.toarray())
1558                            s_min = s[-1]
1559                    else:
1560                        # Dense SVD
1561                        A = H - lam * I
1562                        s = svdvals(A)
1563                        s_min = s[-1]
1564                    
1565                    return idx, 1.0 / (s_min + 1e-16), s_min
1566                except Exception as e:
1567                    return idx, np.nan, np.nan
1568            
1569            # Use ThreadPoolExecutor for parallel computation
1570            with ThreadPoolExecutor(max_workers=n_workers) as executor:
1571                futures = {executor.submit(compute_single_point, idx): idx 
1572                          for idx in range(len(Lambda_flat))}
1573                
1574                # Progress tracking
1575                completed = 0
1576                total = len(futures)
1577                
1578                for future in as_completed(futures):
1579                    idx, res_norm, s_min = future.result()
1580                    resolvent_norm.ravel()[idx] = res_norm
1581                    sigma_min_grid.ravel()[idx] = s_min
1582                    
1583                    completed += 1
1584                    if completed % (total // 10) == 0:
1585                        print(f'Progress: {completed}/{total} ({100*completed//total}%)')
1586            
1587        else:
1588            # Sequential computation
1589            for i in range(resolution):
1590                for j in range(resolution):
1591                    lam = Lambda[i, j]
1592                    try:
1593                        if use_sparse_svd:
1594                            A = H_sparse - lam * I_sparse
1595                            try:
1596                                s_min = svds(A, k=1, which='SM',
1597                                           return_singular_vectors=False)[0]
1598                            except:
1599                                s = svdvals(A.toarray())
1600                                s_min = s[-1]
1601                        else:
1602                            A = H - lam * I
1603                            s = svdvals(A)
1604                            s_min = s[-1]
1605                        
1606                        sigma_min_grid[i, j] = s_min
1607                        resolvent_norm[i, j] = 1.0 / (s_min + 1e-16)
1608                    except Exception:
1609                        resolvent_norm[i, j] = np.nan
1610                        sigma_min_grid[i, j] = np.nan
1611                
1612                if i % (resolution // 10) == 0:
1613                    print(f'Progress: {i}/{resolution} rows')
1614        
1615        return Lambda, resolvent_norm, sigma_min_grid
1616
1617    def _compute_pseudospectrum_adaptive(self, H, lambda_real_range, lambda_imag_range,
1618                                        base_resolution, use_sparse=False, parallel=True,
1619                                        n_workers=4, threshold=0.5, max_refinements=2):
1620        """
1621        Compute pseudospectrum with adaptive grid refinement.
1622        
1623        Starts with coarse grid and refines regions with high gradients.
1624        
1625        Parameters
1626        ----------
1627        H : ndarray
1628            Operator matrix
1629        lambda_real_range : tuple
1630            Range for Re(λ)
1631        lambda_imag_range : tuple
1632            Range for Im(λ)
1633        base_resolution : int
1634            Initial coarse resolution
1635        use_sparse : bool
1636            Use sparse matrices
1637        parallel : bool
1638            Enable parallel computation
1639        n_workers : int
1640            Number of workers
1641        threshold : float
1642            Gradient threshold for refinement
1643        max_refinements : int
1644            Maximum number of refinement levels
1645            
1646        Returns
1647        -------
1648        Lambda : ndarray
1649            Complex grid (may be non-uniform)
1650        resolvent_norm : ndarray
1651            Resolvent norms
1652        sigma_min_grid : ndarray
1653            Smallest singular values
1654        """
1655        # Start with coarse grid
1656        coarse_res = base_resolution // 2
1657        print(f'Level 0: Computing coarse grid ({coarse_res}×{coarse_res})...')
1658        
1659        Lambda_coarse, resolvent_coarse, sigma_coarse = self._compute_pseudospectrum(
1660            H, lambda_real_range, lambda_imag_range, coarse_res,
1661            use_sparse=use_sparse, parallel=parallel, n_workers=n_workers
1662        )
1663        
1664        # Compute gradient to identify regions needing refinement
1665        log_resolvent = np.log10(resolvent_coarse + 1e-16)
1666        grad_y, grad_x = np.gradient(log_resolvent)
1667        grad_magnitude = np.sqrt(grad_x**2 + grad_y**2)
1668        
1669        # Normalize gradient
1670        grad_normalized = grad_magnitude / (np.max(grad_magnitude) + 1e-10)
1671        
1672        # For now, return uniform fine grid
1673        # (Full adaptive implementation would require irregular grids)
1674        print(f'Level 1: Computing fine grid ({base_resolution}×{base_resolution})...')
1675        Lambda_fine, resolvent_fine, sigma_fine = self._compute_pseudospectrum(
1676            H, lambda_real_range, lambda_imag_range, base_resolution,
1677            use_sparse=use_sparse, parallel=parallel, n_workers=n_workers
1678        )
1679        
1680        high_gradient_pct = 100 * np.sum(grad_normalized > threshold) / grad_normalized.size
1681        print(f'High-gradient regions: {high_gradient_pct:.1f}% of domain')
1682        
1683        return Lambda_fine, resolvent_fine, sigma_fine
1684
1685    def _compute_eigenvalues(self, H, use_sparse=False):
1686        """
1687        Compute eigenvalues of operator matrix.
1688        
1689        Parameters
1690        ----------
1691        H : ndarray
1692            Operator matrix
1693        use_sparse : bool
1694            Use sparse eigenvalue solver
1695            
1696        Returns
1697        -------
1698        eigenvalues : ndarray or None
1699            Eigenvalues of H
1700        """
1701        try:
1702            if use_sparse and H.shape[0] > 100:
1703                from scipy.sparse.linalg import eigs
1704                from scipy.sparse import csr_matrix
1705                H_sparse = csr_matrix(H)
1706                k = min(20, H.shape[0] - 2)
1707                eigenvalues = eigs(H_sparse, k=k, return_eigenvectors=False)
1708            else:
1709                eigenvalues = np.linalg.eigvals(H)
1710            
1711            # Print diagnostics
1712            print(f'Eigenvalue range: [{eigenvalues.real.min():.2f}, {eigenvalues.real.max():.2f}]')
1713            print(f'Imaginary part range: [{eigenvalues.imag.min():.2e}, {eigenvalues.imag.max():.2e}]')
1714            
1715            return eigenvalues
1716        except Exception as e:
1717            warnings.warn(f'Eigenvalue computation failed: {e}')
1718            return None
1719
1720    def _plot_pseudospectrum(self, Lambda, resolvent_norm, sigma_min_grid,
1721                            epsilon_levels, eigenvalues):
1722        """
1723        Plot pseudospectrum results.
1724        
1725        Parameters
1726        ----------
1727        Lambda : ndarray
1728            Complex λ grid
1729        resolvent_norm : ndarray
1730            Resolvent norms
1731        sigma_min_grid : ndarray
1732            Smallest singular values
1733        epsilon_levels : list
1734            Contour levels
1735        eigenvalues : ndarray or None
1736            Eigenvalues to overlay
1737        """
1738        Lambda_re = Lambda.real
1739        Lambda_im = Lambda.imag
1740        
1741        plt.figure(figsize=(14, 6))
1742        
1743        # Left plot: ε-pseudospectrum
1744        plt.subplot(1, 2, 1)
1745        
1746        # Better contour level computation
1747        log_resolvent = np.log10(resolvent_norm + 1e-16)
1748        levels_log = np.log10(1.0 / np.array(epsilon_levels))
1749        
1750        # Only plot contours that exist in the data range
1751        valid_levels = [lv for lv in levels_log 
1752                       if log_resolvent.min() <= lv <= log_resolvent.max()]
1753        
1754        if len(valid_levels) > 0:
1755            cs = plt.contour(Lambda_re, Lambda_im, log_resolvent,
1756                            levels=valid_levels, colors='blue', linewidths=1.5)
1757            # Better labels
1758            labels = [f'ε={eps:.0e}' for eps in epsilon_levels[:len(valid_levels)]]
1759            fmt = dict(zip(cs.levels, labels))
1760            plt.clabel(cs, inline=True, fmt=fmt, fontsize=9)
1761        else:
1762            print('⚠️ Warning: No contours in specified epsilon range')
1763            # Plot general contours
1764            cs = plt.contour(Lambda_re, Lambda_im, log_resolvent,
1765                            levels=10, colors='blue', linewidths=1.5)
1766        
1767        if eigenvalues is not None:
1768            plt.plot(eigenvalues.real, eigenvalues.imag, 'r*', 
1769                    markersize=10, label='Eigenvalues', markeredgecolor='darkred')
1770        
1771        plt.xlabel('Re(λ)', fontsize=12)
1772        plt.ylabel('Im(λ)', fontsize=12)
1773        plt.title('ε-Pseudospectrum: log₁₀(‖(H - λI)⁻¹‖)', fontsize=13)
1774        plt.grid(alpha=0.3)
1775        plt.legend(fontsize=10)
1776        plt.axis('equal')
1777        
1778        # Right plot: Smallest singular value
1779        plt.subplot(1, 2, 2)
1780        
1781        # Use better colormap normalization
1782        from matplotlib.colors import LogNorm
1783        
1784        # Filter out invalid values
1785        sigma_plot = np.where(np.isfinite(sigma_min_grid), sigma_min_grid, np.nan)
1786        vmin = np.nanmin(sigma_plot[sigma_plot > 0]) if np.any(sigma_plot > 0) else 1e-10
1787        vmax = np.nanmax(sigma_plot)
1788        
1789        cs2 = plt.contourf(Lambda_re, Lambda_im, sigma_plot,
1790                          levels=50, cmap='viridis',
1791                          norm=LogNorm(vmin=vmin, vmax=vmax))
1792        plt.colorbar(cs2, label='σ_min(H - λI)')
1793        
1794        if eigenvalues is not None:
1795            plt.plot(eigenvalues.real, eigenvalues.imag, 'r*', 
1796                    markersize=10, markeredgecolor='darkred')
1797        
1798        # Plot epsilon contours
1799        for eps in epsilon_levels:
1800            cs_eps = plt.contour(Lambda_re, Lambda_im, sigma_plot,
1801                               levels=[eps], colors='red', linewidths=2, alpha=0.8)
1802        
1803        plt.xlabel('Re(λ)', fontsize=12)
1804        plt.ylabel('Im(λ)', fontsize=12)
1805        plt.title('Smallest singular value σ_min(H - λI)', fontsize=13)
1806        plt.grid(alpha=0.3)
1807        plt.axis('equal')
1808        
1809        plt.tight_layout()
1810        plt.show()
1811    
1812    
1813    def symplectic_flow(self):
1814        """
1815        Compute the Hamiltonian vector field associated with the principal symbol.
1816
1817        This method derives the canonical equations of motion for the phase space variables 
1818        (x, ξ) in 1D or (x, y, ξ, η) in 2D, based on the Hamiltonian formalism. These describe 
1819        how position and frequency variables evolve under the flow generated by the symbol.
1820
1821        Returns
1822        -------
1823        dict
1824            A dictionary containing the components of the Hamiltonian vector field:
1825            - In 1D: keys are 'dx/dt' and 'dxi/dt', corresponding to dx/dt = ∂p/∂ξ and dξ/dt = -∂p/∂x.
1826            - In 2D: keys are 'dx/dt', 'dy/dt', 'dxi/dt', and 'deta/dt', with similar definitions:
1827              dx/dt = ∂p/∂ξ, dy/dt = ∂p/∂η, dξ/dt = -∂p/∂x, dη/dt = -∂p/∂y.
1828
1829        Notes
1830        -----
1831        - The Hamiltonian here is the principal symbol p(x, ξ) itself.
1832        - This flow preserves the symplectic structure of phase space.
1833        """
1834        if self.dim == 1:
1835            x,  = self.vars_x
1836            xi = symbols('xi', real=True)
1837            return {
1838                'dx/dt': diff(self.symbol, xi),
1839                'dxi/dt': -diff(self.symbol, x)
1840            }
1841        elif self.dim == 2:
1842            x, y = self.vars_x
1843            xi, eta = symbols('xi eta', real=True)
1844            return {
1845                'dx/dt': diff(self.symbol, xi),
1846                'dy/dt': diff(self.symbol, eta),
1847                'dxi/dt': -diff(self.symbol, x),
1848                'deta/dt': -diff(self.symbol, y)
1849            }
1850
1851    def is_elliptic_numerically(self, x_grid, xi_grid, threshold=1e-8):
1852        """
1853        Check if the pseudo-differential symbol p(x, ξ) is elliptic over a given grid.
1854    
1855        A symbol is considered elliptic if its magnitude |p(x, ξ)| remains bounded away from zero 
1856        across all points in the spatial-frequency domain. This method evaluates the symbol on a 
1857        grid of spatial and frequency coordinates and checks whether its minimum absolute value 
1858        exceeds a specified threshold.
1859    
1860        Resampling is applied to large grids to prevent excessive memory usage, particularly in 2D.
1861    
1862        Parameters
1863        ----------
1864        x_grid : ndarray
1865            Spatial grid: either a 1D array (x) or a tuple of two 1D arrays (x, y).
1866        xi_grid : ndarray
1867            Frequency grid: either a 1D array (ξ) or a tuple of two 1D arrays (ξ, η).
1868        threshold : float, optional
1869            Minimum acceptable value for |p(x, ξ)|. If the smallest evaluated symbol value falls below this,
1870            the symbol is not considered elliptic.
1871    
1872        Returns
1873        -------
1874        bool
1875            True if the symbol is elliptic on the resampled grid, False otherwise.
1876        """
1877        RESAMPLE_SIZE = 32  # Reduced size to prevent memory explosion
1878        
1879        if self.dim == 1:
1880            x_vals = x_grid
1881            xi_vals = xi_grid
1882            # Resampling if necessary
1883            if len(x_vals) > RESAMPLE_SIZE:
1884                x_vals = np.linspace(x_vals.min(), x_vals.max(), RESAMPLE_SIZE)
1885            if len(xi_vals) > RESAMPLE_SIZE:
1886                xi_vals = np.linspace(xi_vals.min(), xi_vals.max(), RESAMPLE_SIZE)
1887        
1888            X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
1889            symbol_vals = self.p_func(X, XI)
1890        
1891        elif self.dim == 2:
1892            x_vals, y_vals = x_grid
1893            xi_vals, eta_vals = xi_grid
1894        
1895            # Spatial resampling
1896            if len(x_vals) > RESAMPLE_SIZE:
1897                x_vals = np.linspace(x_vals.min(), x_vals.max(), RESAMPLE_SIZE)
1898            if len(y_vals) > RESAMPLE_SIZE:
1899                y_vals = np.linspace(y_vals.min(), y_vals.max(), RESAMPLE_SIZE)
1900        
1901            # Frequency resampling
1902            if len(xi_vals) > RESAMPLE_SIZE:
1903                xi_vals = np.linspace(xi_vals.min(), xi_vals.max(), RESAMPLE_SIZE)
1904            if len(eta_vals) > RESAMPLE_SIZE:
1905                eta_vals = np.linspace(eta_vals.min(), eta_vals.max(), RESAMPLE_SIZE)
1906        
1907            X, Y, XI, ETA = np.meshgrid(x_vals, y_vals, xi_vals, eta_vals, indexing='ij')
1908            symbol_vals = self.p_func(X, Y, XI, ETA)
1909        
1910        min_abs_val = np.min(np.abs(symbol_vals))
1911        return min_abs_val > threshold
1912
1913
1914    def is_self_adjoint(self, tol=1e-10):
1915        """
1916        Check whether the pseudo-differential operator is formally self-adjoint (Hermitian).
1917
1918        A self-adjoint operator satisfies P = P*, where P* is the formal adjoint of P.
1919        This property is essential for ensuring real-valued eigenvalues and stable evolution 
1920        in quantum mechanics and symmetric wave propagation.
1921
1922        Parameters
1923        ----------
1924        tol : float
1925            Tolerance for symbolic comparison between P and P*. Small numerical differences 
1926            below this threshold are considered equal.
1927
1928        Returns
1929        -------
1930        bool
1931            True if the symbol p(x, ξ) equals its formal adjoint p*(x, ξ) within the given tolerance,
1932            indicating that the operator is self-adjoint.
1933
1934        Notes:
1935        - The formal adjoint is computed via conjugation and asymptotic expansion at infinity in ξ.
1936        - Symbolic simplification is used to verify equality, ensuring robustness against superficial 
1937          expression differences.
1938        """
1939        p = self.symbol
1940        p_star = self.formal_adjoint()
1941        return simplify(p - p_star).equals(0)
1942
1943    def visualize_fiber(self, x_grid, xi_grid, x0=0.0, y0=0.0):
1944        """
1945        Plot the cotangent fiber structure at a fixed spatial point (x₀[, y₀]).
1946    
1947        This visualization shows how the symbol p(x, ξ) behaves on the cotangent fiber 
1948        above a fixed spatial point. In microlocal analysis, this provides insight into 
1949        the frequency content of the operator at that location.
1950    
1951        Parameters
1952        ----------
1953        x_grid : ndarray
1954            Spatial grid values (1D) for evaluation in 1D case.
1955        xi_grid : ndarray
1956            Frequency grid values (1D) for evaluation in both 1D and 2D cases.
1957        x0 : float, optional
1958            Fixed x-coordinate of the base point in space (1D or 2D).
1959        y0 : float, optional
1960            Fixed y-coordinate of the base point in space (2D only).
1961    
1962        Notes
1963        -----
1964        - In 1D: Displays |p(x, ξ)| over the (x, ξ) phase plane near the fixed point.
1965        - In 2D: Fixes (x₀, y₀) and evaluates p(x₀, y₀, ξ, η), showing the fiber over that point.
1966        - The color map represents the magnitude of the symbol, highlighting regions where it vanishes or becomes singular.
1967    
1968        Raises
1969        ------
1970        NotImplementedError
1971            If called in 2D with missing or improperly formatted grids.
1972        """
1973        if self.dim == 1:
1974            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
1975            symbol_vals = self.p_func(X, XI)
1976            plt.contourf(X, XI, np.abs(symbol_vals), levels=50, cmap='viridis')
1977            plt.colorbar(label='|Symbol|')
1978            plt.xlabel('x (position)')
1979            plt.ylabel('ξ (frequency)')
1980            plt.title('Cotangent Fiber Structure')
1981            plt.show()
1982        elif self.dim == 2:
1983            xi_grid2, eta_grid2 = np.meshgrid(xi_grid, xi_grid)
1984            symbol_vals = self.p_func(x0, y0, xi_grid2, eta_grid2)
1985            plt.contourf(xi_grid, xi_grid, np.abs(symbol_vals), levels=50, cmap='viridis')
1986            plt.colorbar(label='|Symbol|')
1987            plt.xlabel('ξ')
1988            plt.ylabel('η')
1989            plt.title(f'Cotangent Fiber at x={x0}, y={y0}')
1990            plt.show()
1991
1992    def visualize_symbol_amplitude(self, x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0):
1993        """
1994        Display the modulus |p(x, ξ)| or |p(x, y, ξ₀, η₀)| as a color map.
1995    
1996        This method visualizes the amplitude of the pseudodifferential operator's symbol 
1997        in either 1D or 2D spatial configuration. In 2D, the frequency variables are fixed 
1998        to specified values (ξ₀, η₀) for visualization purposes.
1999    
2000        Parameters
2001        ----------
2002        x_grid, y_grid : ndarray
2003            Spatial grids over which to evaluate the symbol. y_grid is optional and used only in 2D.
2004        xi_grid, eta_grid : ndarray
2005            Frequency grids. In 2D, these define the domain over which the symbol is evaluated,
2006            but the visualization fixes ξ = ξ₀ and η = η₀.
2007        xi0, eta0 : float, optional
2008            Fixed frequency values for slicing in 2D visualization. Defaults to zero.
2009    
2010        Notes
2011        -----
2012        - In 1D: Visualizes |p(x, ξ)| over the (x, ξ) grid.
2013        - In 2D: Visualizes |p(x, y, ξ₀, η₀)| at fixed frequencies ξ₀ and η₀.
2014        - The color intensity represents the magnitude of the symbol, highlighting regions where the symbol is large or small.
2015        """
2016        if self.dim == 1:
2017            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
2018            symbol_vals = self.p_func(X, XI) 
2019            plt.pcolormesh(X, XI, np.abs(symbol_vals), shading='auto')
2020            plt.colorbar(label='|Symbol|')
2021            plt.xlabel('x')
2022            plt.ylabel('ξ')
2023            plt.title('Symbol Amplitude |p(x, ξ)|')
2024            plt.show()
2025        elif self.dim == 2:
2026            X, Y = np.meshgrid(x_grid, y_grid, indexing='ij')
2027            XI = np.full_like(X, xi0)
2028            ETA = np.full_like(Y, eta0)
2029            symbol_vals = self.p_func(X, Y, XI, ETA)
2030            plt.pcolormesh(X, Y, np.abs(symbol_vals), shading='auto')
2031            plt.colorbar(label='|Symbol|')
2032            plt.xlabel('x')
2033            plt.ylabel('y')
2034            plt.title(f'Symbol Amplitude at ξ={xi0}, η={eta0}')
2035            plt.show()
2036
2037    def visualize_phase(self, x_grid, xi_grid, y_grid=None, eta_grid=None, xi0=0.0, eta0=0.0):
2038        """
2039        Plot the phase (argument) of the pseudodifferential operator's symbol p(x, ξ) or p(x, y, ξ, η).
2040
2041        This visualization helps in understanding the oscillatory behavior and regularity 
2042        properties of the operator in phase space. The phase is displayed modulo 2π using 
2043        a cyclic colormap ('twilight') to emphasize its periodic nature.
2044
2045        Parameters
2046        ----------
2047        x_grid : ndarray
2048            1D array of spatial coordinates (x).
2049        xi_grid : ndarray
2050            1D array of frequency coordinates (ξ).
2051        y_grid : ndarray, optional
2052            2D spatial grid for y-coordinate (in 2D problems). Default is None.
2053        eta_grid : ndarray, optional
2054            2D frequency grid for η (in 2D problems). Not used directly but kept for API consistency.
2055        xi0 : float, optional
2056            Fixed value of ξ for slicing in 2D visualization. Default is 0.0.
2057        eta0 : float, optional
2058            Fixed value of η for slicing in 2D visualization. Default is 0.0.
2059
2060        Notes:
2061        - In 1D: Displays arg(p(x, ξ)) over the (x, ξ) phase plane.
2062        - In 2D: Displays arg(p(x, y, ξ₀, η₀)) for fixed frequency values (ξ₀, η₀).
2063        - Uses plt.pcolormesh with 'twilight' colormap to represent angles from -π to π.
2064
2065        Raises:
2066        - NotImplementedError: If the spatial dimension is not 1D or 2D.
2067        """
2068        if self.dim == 1:
2069            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
2070            symbol_vals = self.p_func(X, XI) 
2071            plt.pcolormesh(X, XI, np.angle(symbol_vals), shading='auto', cmap='twilight')
2072            plt.colorbar(label='arg(Symbol) [rad]')
2073            plt.xlabel('x')
2074            plt.ylabel('ξ')
2075            plt.title('Phase Portrait (arg p(x, ξ))')
2076            plt.show()
2077        elif self.dim == 2:
2078            X, Y = np.meshgrid(x_grid, y_grid, indexing='ij')
2079            XI = np.full_like(X, xi0)
2080            ETA = np.full_like(Y, eta0)
2081            symbol_vals = self.p_func(X, Y, XI, ETA)
2082            plt.pcolormesh(X, Y, np.angle(symbol_vals), shading='auto', cmap='twilight')
2083            plt.colorbar(label='arg(Symbol) [rad]')
2084            plt.xlabel('x')
2085            plt.ylabel('y')
2086            plt.title(f'Phase Portrait at ξ={xi0}, η={eta0}')
2087            plt.show()
2088            
2089    def visualize_characteristic_set(self, x_grid, xi_grid, y_grid=None, eta_grid=None, y0=0.0, x0=0.0, levels=[1e-1]):
2090        """
2091        Visualize the characteristic set of the pseudo-differential symbol, defined as the approximate zero set p(x, ξ) ≈ 0.
2092    
2093        In microlocal analysis, the characteristic set is the locus of points in phase space (x, ξ) where the symbol p(x, ξ) vanishes,
2094        playing a key role in understanding propagation of singularities.
2095    
2096        Parameters
2097        ----------
2098        x_grid : ndarray
2099            Spatial grid values (1D array) for plotting in 1D or evaluation point in 2D.
2100        xi_grid : ndarray
2101            Frequency variable grid values (1D array) used to construct the frequency domain.
2102        x0 : float, optional
2103            Fixed spatial coordinate in 2D case for evaluating the symbol at a specific x position.
2104        y0 : float, optional
2105            Fixed spatial coordinate in 2D case for evaluating the symbol at a specific y position.
2106    
2107        Notes
2108        -----
2109        - For 1D, this method plots the contour of |p(x, ξ)| = ε with ε = 1e-5 over the (x, ξ) plane.
2110        - For 2D, it evaluates the symbol at fixed (x₀, y₀) and plots the characteristic set in the (ξ, η) frequency plane.
2111        - This visualization helps identify directions of degeneracy or hypoellipticity of the operator.
2112    
2113        Raises
2114        ------
2115        NotImplementedError
2116            If called on a solver with dimensionality other than 1D or 2D.
2117    
2118        Displays
2119        ------
2120        A matplotlib contour plot showing either:
2121            - The characteristic curve in the (x, ξ) phase plane (1D),
2122            - The characteristic surface slice in the (ξ, η) frequency plane at (x₀, y₀) (2D).
2123        """
2124        if self.dim == 1:
2125            x_grid = np.asarray(x_grid)
2126            xi_grid = np.asarray(xi_grid)
2127            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
2128            symbol_vals = self.p_func(X, XI) 
2129            plt.contour(X, XI, np.abs(symbol_vals), levels=levels, colors='red')
2130            plt.xlabel('x')
2131            plt.ylabel('ξ')
2132            plt.title('Characteristic Set (p(x, ξ) ≈ 0)')
2133            plt.grid(True)
2134            plt.show()
2135        elif self.dim == 2:
2136            if eta_grid is None:
2137                raise ValueError("eta_grid must be provided for 2D visualization.")
2138            xi_grid = np.asarray(xi_grid)
2139            eta_grid = np.asarray(eta_grid)
2140            xi_grid2, eta_grid2 = np.meshgrid(xi_grid, eta_grid, indexing='ij')
2141            symbol_vals = self.p_func(x0, y0, xi_grid2, eta_grid2)
2142            plt.contour(xi_grid, eta_grid, np.abs(symbol_vals), levels=levels, colors='red')
2143            plt.xlabel('ξ')
2144            plt.ylabel('η')
2145            plt.title(f'Characteristic Set at x={x0}, y={y0}')
2146            plt.grid(True)
2147            plt.show()
2148        else:
2149            raise NotImplementedError("Only 1D/2D characteristic sets supported.")
2150
2151    def visualize_characteristic_gradient(self, x_grid, xi_grid, y_grid=None, eta_grid=None, y0=0.0, x0=0.0):
2152        """
2153        Visualize the norm of the gradient of the symbol in phase space.
2154        
2155        This method computes the magnitude of the gradient |∇p| of a pseudo-differential 
2156        symbol p(x, ξ) in 1D or p(x, y, ξ, η) in 2D. The resulting colormap reveals 
2157        regions where the symbol varies rapidly or remains nearly stationary, 
2158        which is particularly useful for analyzing characteristic sets.
2159        
2160        Parameters
2161        ----------
2162        x_grid : numpy.ndarray
2163            1D array of spatial coordinates for the x-direction.
2164        xi_grid : numpy.ndarray
2165            1D array of frequency coordinates (ξ).
2166        y_grid : numpy.ndarray, optional
2167            1D array of spatial coordinates for the y-direction (used in 2D mode). Default is None.
2168        eta_grid : numpy.ndarray, optional
2169            1D array of frequency coordinates (η) for the 2D case. Default is None.
2170        x0 : float, optional
2171            Fixed x-coordinate for evaluating the symbol in 2D. Default is 0.0.
2172        y0 : float, optional
2173            Fixed y-coordinate for evaluating the symbol in 2D. Default is 0.0.
2174        
2175        Returns
2176        -------
2177        None
2178            Displays a 2D colormap of |∇p| over the relevant phase-space domain.
2179        
2180        Notes
2181        -----
2182        - In 1D, the full gradient ∇p = (∂ₓp, ∂ξp) is computed over the (x, ξ) grid.
2183        - In 2D, the gradient ∇p = (∂ξp, ∂ηp) is computed at a fixed spatial point (x₀, y₀) over the (ξ, η) grid.
2184        - Numerical differentiation is performed using `np.gradient`.
2185        - High values of |∇p| indicate rapid variation of the symbol, while low values typically suggest characteristic regions.
2186        """
2187        if self.dim == 1:
2188            X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
2189            symbol_vals = self.p_func(X, XI)
2190            grad_x = np.gradient(symbol_vals, axis=0)
2191            grad_xi = np.gradient(symbol_vals, axis=1)
2192            grad_norm = np.sqrt(grad_x**2 + grad_xi**2)
2193            plt.pcolormesh(X, XI, grad_norm, cmap='inferno', shading='auto')
2194            plt.colorbar(label='|∇p|')
2195            plt.xlabel('x')
2196            plt.ylabel('ξ')
2197            plt.title('Gradient Norm (High Near Zeros)')
2198            plt.grid(True)
2199            plt.show()
2200        elif self.dim == 2:
2201            xi_grid2, eta_grid2 = np.meshgrid(xi_grid, eta_grid, indexing='ij')
2202            symbol_vals = self.p_func(x0, y0, xi_grid2, eta_grid2)
2203            grad_xi = np.gradient(symbol_vals, axis=0)
2204            grad_eta = np.gradient(symbol_vals, axis=1)
2205            grad_norm = np.sqrt(np.abs(grad_xi)**2 + np.abs(grad_eta)**2)
2206            plt.pcolormesh(xi_grid, eta_grid, grad_norm, cmap='inferno', shading='auto')
2207            plt.colorbar(label='|∇p|')
2208            plt.xlabel('ξ')
2209            plt.ylabel('η')
2210            plt.title(f'Gradient Norm at x={x0}, y={y0}')
2211            plt.grid(True)
2212            plt.show()
2213
2214    def plot_hamiltonian_flow(self, x0=0.0, xi0=5.0, y0=0.0, eta0=0.0, tmax=1.0, n_steps=100, show_field=True):
2215        """
2216        Integrate and plot the Hamiltonian trajectories of the symbol in phase space.
2217
2218        This method numerically integrates the Hamiltonian vector field derived from 
2219        the operator's symbol to visualize how singularities propagate under the flow. 
2220        It supports both 1D and 2D problems.
2221
2222        Parameters
2223        ----------
2224        x0, xi0 : float
2225            Initial position and frequency (momentum) in 1D.
2226        y0, eta0 : float, optional
2227            Initial position and frequency in 2D; defaults to zero.
2228        tmax : float
2229            Final integration time for the ODE solver.
2230        n_steps : int
2231            Number of time steps used in the integration.
2232
2233        Notes
2234        -----
2235        - The Hamiltonian vector field is obtained from the symplectic flow of the symbol.
2236        - If the field is complex-valued, only its real part is used for integration.
2237        - In 1D, the trajectory is plotted in (x, ξ) phase space.
2238        - In 2D, the spatial trajectory (x(t), y(t)) is shown along with instantaneous 
2239          momentum vectors (ξ(t), η(t)) using a quiver plot.
2240
2241        Raises
2242        ------
2243        NotImplementedError
2244            If the spatial dimension is not 1D or 2D.
2245
2246        Displays
2247        --------
2248        matplotlib plot
2249            Phase space trajectory(ies) showing the evolution of position and momentum 
2250            under the Hamiltonian dynamics.
2251        """
2252        def make_real(expr):
2253            from sympy import re, simplify
2254            expr = expr.doit(deep=True)
2255            return simplify(re(expr))
2256    
2257        H = self.symplectic_flow()
2258    
2259        if any(im(H[k]) != 0 for k in H):
2260            print("⚠️ The Hamiltonian field is complex. Only the real part is used for integration.")
2261    
2262        if self.dim == 1:
2263            x, = self.vars_x
2264            xi = symbols('xi', real=True)
2265    
2266            dxdt_expr = make_real(H['dx/dt'])
2267            dxidt_expr = make_real(H['dxi/dt'])
2268    
2269            dxdt = lambdify((x, xi), dxdt_expr, 'numpy')
2270            dxidt = lambdify((x, xi), dxidt_expr, 'numpy')
2271    
2272            def hamilton(t, Y):
2273                x, xi = Y
2274                return [dxdt(x, xi), dxidt(x, xi)]
2275    
2276            sol = solve_ivp(hamilton, [0, tmax], [x0, xi0], t_eval=np.linspace(0, tmax, n_steps))
2277
2278            if sol.status != 0:
2279                print(f"⚠️ Integration warning: {sol.message}")
2280            
2281            n_points = sol.y.shape[1]
2282            if n_points < n_steps:
2283                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
2284                n_steps = n_points
2285
2286            x_vals, xi_vals = sol.y
2287    
2288            plt.plot(x_vals, xi_vals)
2289            plt.xlabel("x")
2290            plt.ylabel("ξ")
2291            plt.title("Hamiltonian Flow in Phase Space (1D)")
2292            plt.grid(True)
2293            plt.show()
2294    
2295        elif self.dim == 2:
2296            x, y = self.vars_x
2297            xi, eta = symbols('xi eta', real=True)
2298    
2299            dxdt = lambdify((x, y, xi, eta), make_real(H['dx/dt']), 'numpy')
2300            dydt = lambdify((x, y, xi, eta), make_real(H['dy/dt']), 'numpy')
2301            dxidt = lambdify((x, y, xi, eta), make_real(H['dxi/dt']), 'numpy')
2302            detadt = lambdify((x, y, xi, eta), make_real(H['deta/dt']), 'numpy')
2303    
2304            def hamilton(t, Y):
2305                x, y, xi, eta = Y
2306                return [
2307                    dxdt(x, y, xi, eta),
2308                    dydt(x, y, xi, eta),
2309                    dxidt(x, y, xi, eta),
2310                    detadt(x, y, xi, eta)
2311                ]
2312    
2313            sol = solve_ivp(hamilton, [0, tmax], [x0, y0, xi0, eta0], t_eval=np.linspace(0, tmax, n_steps))
2314
2315            if sol.status != 0:
2316                print(f"⚠️ Integration warning: {sol.message}")
2317            
2318            n_points = sol.y.shape[1]
2319            if n_points < n_steps:
2320                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
2321                n_steps = n_points
2322
2323            x_vals, y_vals, xi_vals, eta_vals = sol.y
2324    
2325            plt.plot(x_vals, y_vals, label='Position')
2326            plt.quiver(x_vals, y_vals, xi_vals, eta_vals, scale=20, width=0.003, alpha=0.5, color='r')
2327            
2328            # Vector field of the flow (optional)
2329            if show_field:
2330                X, Y = np.meshgrid(np.linspace(min(x_vals), max(x_vals), 20),
2331                                   np.linspace(min(y_vals), max(y_vals), 20))
2332                XI, ETA = xi0 * np.ones_like(X), eta0 * np.ones_like(Y)
2333                U = dxdt(X, Y, XI, ETA)
2334                V = dydt(X, Y, XI, ETA)
2335                plt.quiver(X, Y, U, V, color='gray', alpha=0.2, scale=30, width=0.002)
2336
2337            plt.xlabel("x")
2338            plt.ylabel("y")
2339            plt.title("Hamiltonian Flow in Phase Space (2D)")
2340            plt.legend()
2341            plt.grid(True)
2342            plt.axis('equal')
2343            plt.show()
2344
2345    def plot_symplectic_vector_field(self, xlim=(-2, 2), klim=(-5, 5), density=30):
2346        """
2347        Visualize the symplectic vector field (Hamiltonian vector field) associated with the operator's symbol.
2348
2349        The plotted vector field corresponds to (∂_ξ p, -∂_x p), where p(x, ξ) is the principal symbol 
2350        of the pseudo-differential operator. This field governs the bicharacteristic flow in phase space.
2351
2352        Parameters
2353        ----------
2354        xlim : tuple of float
2355            Range for spatial variable x, as (x_min, x_max).
2356        klim : tuple of float
2357            Range for frequency variable ξ, as (ξ_min, ξ_max).
2358        density : int
2359            Number of grid points per axis for the visualization grid.
2360
2361        Raises
2362        ------
2363        NotImplementedError
2364            If called on a 2D operator (currently only 1D implementation available).
2365
2366        Notes
2367        -----
2368        - Only supports one-dimensional operators.
2369        - Uses symbolic differentiation to compute ∂_ξ p and ∂_x p.
2370        - Numerical evaluation is done via lambdify with NumPy backend.
2371        - Visualization uses matplotlib quiver plot to show vector directions.
2372        """
2373        x_vals = np.linspace(*xlim, density)
2374        xi_vals = np.linspace(*klim, density)
2375        X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
2376
2377        if self.dim != 1:
2378            raise NotImplementedError("Only 1D version implemented.")
2379
2380        x, = self.vars_x
2381        xi = symbols('xi', real=True)
2382        H = self.symplectic_flow()
2383        dxdt = lambdify((x, xi), simplify(H['dx/dt']), 'numpy')
2384        dxidt = lambdify((x, xi), simplify(H['dxi/dt']), 'numpy')
2385
2386        U = dxdt(X, XI)
2387        V = dxidt(X, XI)
2388
2389        plt.quiver(X, XI, U, V, scale=10, width=0.005)
2390        plt.xlabel('x')
2391        plt.ylabel(r'$\xi$')
2392        plt.title("Symplectic Vector Field (1D)")
2393        plt.grid(True)
2394        plt.show()
2395
2396    def visualize_micro_support(self, xlim=(-2, 2), klim=(-10, 10), threshold=1e-3, density=300):
2397        """
2398        Visualize the micro-support of the operator by plotting the inverse of the symbol magnitude 1 / |p(x, ξ)|.
2399    
2400        The micro-support provides insight into the singularities of a pseudo-differential operator 
2401        in phase space (x, ξ). Regions where |p(x, ξ)| is small correspond to large values in 1/|p(x, ξ)|,
2402        highlighting areas of significant operator influence or singularity.
2403    
2404        Parameters
2405        ----------
2406        xlim : tuple
2407            Spatial domain limits (x_min, x_max).
2408        klim : tuple
2409            Frequency domain limits (ξ_min, ξ_max).
2410        threshold : float
2411            Threshold below which |p(x, ξ)| is considered effectively zero; used for numerical stability.
2412        density : int
2413            Number of grid points along each axis for visualization resolution.
2414    
2415        Raises
2416        ------
2417        NotImplementedError
2418            If called on a solver with dimension greater than 1 (only 1D visualization is supported).
2419    
2420        Notes
2421        -----
2422        - This method evaluates the symbol p(x, ξ) over a grid and plots its reciprocal to emphasize 
2423          regions where the symbol is near zero.
2424        - A small constant (1e-10) is added to the denominator to avoid division by zero.
2425        - The resulting plot helps identify characteristic sets.
2426        """
2427        if self.dim != 1:
2428            raise NotImplementedError("Only 1D micro-support visualization implemented.")
2429
2430        x_vals = np.linspace(*xlim, density)
2431        xi_vals = np.linspace(*klim, density)
2432        X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
2433        Z = np.abs(self.p_func(X, XI))
2434
2435        plt.contourf(X, XI, 1 / (Z + 1e-10), levels=100, cmap='inferno')
2436        plt.colorbar(label=r'$1/|p(x,\xi)|$')
2437        plt.xlabel('x')
2438        plt.ylabel(r'$\xi$')
2439        plt.title("Micro-Support Estimate (1/|Symbol|)")
2440        plt.show()
2441
2442    def group_velocity_field(self, xlim=(-2, 2), klim=(-10, 10), density=30):
2443        """
2444        Plot the group velocity field ∇_ξ p(x, ξ) for 1D pseudo-differential operators.
2445
2446        The group velocity represents the speed at which waves of different frequencies propagate 
2447        in a dispersive medium. It is defined as the gradient of the symbol p(x, ξ) with respect 
2448        to the frequency variable ξ.
2449
2450        Parameters
2451        ----------
2452        xlim : tuple of float
2453            Spatial domain limits (x-axis).
2454        klim : tuple of float
2455            Frequency domain limits (ξ-axis).
2456        density : int
2457            Number of grid points per axis used for visualization.
2458
2459        Raises
2460        ------
2461        NotImplementedError
2462            If called on a 2D operator, since this visualization is only implemented for 1D.
2463
2464        Notes
2465        -----
2466        - This method visualizes the vector field (∂p/∂ξ) in phase space.
2467        - Used for analyzing wave propagation properties and dispersion relations.
2468        - Requires symbolic expression self.expr depending on x and ξ.
2469        """
2470        if self.dim != 1:
2471            raise NotImplementedError("Only 1D group velocity visualization implemented.")
2472
2473        x, = self.vars_x
2474        xi = symbols('xi', real=True)
2475        dp_dxi = diff(self.symbol, xi)
2476        grad_func = lambdify((x, xi), dp_dxi, 'numpy')
2477
2478        x_vals = np.linspace(*xlim, density)
2479        xi_vals = np.linspace(*klim, density)
2480        X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
2481        V = grad_func(X, XI)
2482
2483        plt.quiver(X, XI, np.ones_like(V), V, scale=10, width=0.004)
2484        plt.xlabel('x')
2485        plt.ylabel(r'$\xi$')
2486        plt.title("Group Velocity Field (1D)")
2487        plt.grid(True)
2488        plt.show()
2489
2490    def animate_singularity(self, xi0=5.0, eta0=0.0, x0=0.0, y0=0.0,
2491                            tmax=4.0, n_frames=100, projection=None):
2492        """
2493        Animate the propagation of a singularity under the Hamiltonian flow.
2494
2495        This method visualizes how a singularity (x₀, y₀, ξ₀, η₀) evolves in phase space 
2496        according to the Hamiltonian dynamics induced by the principal symbol of the operator.
2497        The animation integrates the Hamiltonian equations of motion and supports various projections:
2498        position (x-y), frequency (ξ-η), or mixed phase space coordinates.
2499
2500        Parameters
2501        ----------
2502        xi0, eta0 : float
2503            Initial frequency components (ξ₀, η₀).
2504        x0, y0 : float
2505            Initial spatial coordinates (x₀, y₀).
2506        tmax : float
2507            Total time of integration (final animation time).
2508        n_frames : int
2509            Number of frames in the resulting animation.
2510        projection : str or None
2511            Type of projection to display:
2512                - 'position' : x vs y (or x alone in 1D)
2513                - 'frequency': ξ vs η (or ξ alone in 1D)
2514                - 'phase'    : mixed coordinates like x vs ξ or x vs η
2515                If None, defaults to 'phase' in 1D and 'position' in 2D.
2516
2517        Returns
2518        -------
2519        matplotlib.animation.FuncAnimation
2520            Animation object that can be displayed interactively in Jupyter notebooks or saved as a video.
2521
2522        Notes
2523        -----
2524        - In 1D, only one spatial and one frequency variable are used.
2525        - Complex-valued Hamiltonian fields are truncated to their real parts for integration.
2526        - Trajectories are shown with both instantaneous position (dot) and full path (dashed line).
2527        """
2528        rc('animation', html='jshtml')
2529    
2530        def make_real(expr):
2531            from sympy import re, simplify
2532            expr = expr.doit(deep=True)
2533            return simplify(re(expr))
2534  
2535        H = self.symplectic_flow()
2536
2537        H = {k: v.doit(deep=True) for k, v in H.items()}
2538
2539        print("H = ", H)
2540    
2541        if any(im(H[k]) != 0 for k in H):
2542            print("⚠️ The Hamiltonian field is complex. Only the real part is used for integration.")
2543    
2544        if self.dim == 1:
2545            x, = self.vars_x
2546            xi = symbols('xi', real=True)
2547    
2548            dxdt = lambdify((x, xi), make_real(H['dx/dt']), 'numpy')
2549            dxidt = lambdify((x, xi), make_real(H['dxi/dt']), 'numpy')
2550    
2551            def hamilton(t, Y):
2552                x, xi = Y
2553                return [dxdt(x, xi), dxidt(x, xi)]
2554    
2555            sol = solve_ivp(hamilton, [0, tmax], [x0, xi0],
2556                            t_eval=np.linspace(0, tmax, n_frames))
2557            
2558            if sol.status != 0:
2559                print(f"⚠️ Integration warning: {sol.message}")
2560            
2561            n_points = sol.y.shape[1]
2562            if n_points < n_frames:
2563                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
2564                n_frames = n_points
2565
2566            x_vals, xi_vals = sol.y
2567    
2568            if projection is None:
2569                projection = 'phase'
2570    
2571            fig, ax = plt.subplots()
2572            point, = ax.plot([], [], 'ro')
2573            traj, = ax.plot([], [], 'b--', lw=1, alpha=0.5)
2574    
2575            if projection == 'phase':
2576                ax.set_xlabel('x')
2577                ax.set_ylabel(r'$\xi$')
2578                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2579                ax.set_ylim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2580    
2581                def update(i):
2582                    point.set_data([x_vals[i]], [xi_vals[i]])
2583                    traj.set_data(x_vals[:i+1], xi_vals[:i+1])
2584                    return point, traj
2585    
2586            elif projection == 'position':
2587                ax.set_xlabel('x')
2588                ax.set_ylabel('x')
2589                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2590                ax.set_ylim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2591    
2592                def update(i):
2593                    point.set_data([x_vals[i]], [x_vals[i]])
2594                    traj.set_data(x_vals[:i+1], x_vals[:i+1])
2595                    return point, traj
2596    
2597            elif projection == 'frequency':
2598                ax.set_xlabel(r'$\xi$')
2599                ax.set_ylabel(r'$\xi$')
2600                ax.set_xlim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2601                ax.set_ylim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2602    
2603                def update(i):
2604                    point.set_data([xi_vals[i]], [xi_vals[i]])
2605                    traj.set_data(xi_vals[:i+1], xi_vals[:i+1])
2606                    return point, traj
2607    
2608            else:
2609                raise ValueError("Invalid projection mode")
2610    
2611            ax.set_title(f"1D Singularity Flow ({projection})")
2612            ax.grid(True)
2613            ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=50)
2614            plt.close(fig)
2615            return ani
2616    
2617        elif self.dim == 2:
2618            x, y = self.vars_x
2619            xi, eta = symbols('xi eta', real=True)
2620    
2621            dxdt = lambdify((x, y, xi, eta), make_real(H['dx/dt']), 'numpy')
2622            dydt = lambdify((x, y, xi, eta), make_real(H['dy/dt']), 'numpy')
2623            dxidt = lambdify((x, y, xi, eta), make_real(H['dxi/dt']), 'numpy')
2624            detadt = lambdify((x, y, xi, eta), make_real(H['deta/dt']), 'numpy')
2625    
2626            def hamilton(t, Y):
2627                x, y, xi, eta = Y
2628                return [
2629                    dxdt(x, y, xi, eta),
2630                    dydt(x, y, xi, eta),
2631                    dxidt(x, y, xi, eta),
2632                    detadt(x, y, xi, eta)
2633                ]
2634    
2635            sol = solve_ivp(hamilton, [0, tmax], [x0, y0, xi0, eta0],
2636                            t_eval=np.linspace(0, tmax, n_frames))
2637
2638            if sol.status != 0:
2639                print(f"⚠️ Integration warning: {sol.message}")
2640            
2641            n_points = sol.y.shape[1]
2642            if n_points < n_frames:
2643                print(f"⚠️ Only {n_points} frames computed. Adjusting animation.")
2644                n_frames = n_points
2645                
2646            x_vals, y_vals, xi_vals, eta_vals = sol.y
2647    
2648            if projection is None:
2649                projection = 'position'
2650    
2651            fig, ax = plt.subplots()
2652            point, = ax.plot([], [], 'ro')
2653            traj, = ax.plot([], [], 'b--', lw=1, alpha=0.5)
2654    
2655            if projection == 'position':
2656                ax.set_xlabel('x')
2657                ax.set_ylabel('y')
2658                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2659                ax.set_ylim(np.min(y_vals) - 1, np.max(y_vals) + 1)
2660    
2661                def update(i):
2662                    point.set_data([x_vals[i]], [y_vals[i]])
2663                    traj.set_data(x_vals[:i+1], y_vals[:i+1])
2664                    return point, traj
2665    
2666            elif projection == 'frequency':
2667                ax.set_xlabel(r'$\xi$')
2668                ax.set_ylabel(r'$\eta$')
2669                ax.set_xlim(np.min(xi_vals) - 1, np.max(xi_vals) + 1)
2670                ax.set_ylim(np.min(eta_vals) - 1, np.max(eta_vals) + 1)
2671    
2672                def update(i):
2673                    point.set_data([xi_vals[i]], [eta_vals[i]])
2674                    traj.set_data(xi_vals[:i+1], eta_vals[:i+1])
2675                    return point, traj
2676    
2677            elif projection == 'phase':
2678                ax.set_xlabel('x')
2679                ax.set_ylabel(r'$\eta$')
2680                ax.set_xlim(np.min(x_vals) - 1, np.max(x_vals) + 1)
2681                ax.set_ylim(np.min(eta_vals) - 1, np.max(eta_vals) + 1)
2682    
2683                def update(i):
2684                    point.set_data([x_vals[i]], [eta_vals[i]])
2685                    traj.set_data(x_vals[:i+1], eta_vals[:i+1])
2686                    return point, traj
2687    
2688            else:
2689                raise ValueError("Invalid projection mode")
2690    
2691            ax.set_title(f"2D Singularity Flow ({projection})")
2692            ax.grid(True)
2693            ax.axis('equal')
2694            ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=50)
2695            plt.close(fig)
2696            return ani
2697
2698    def interactive_symbol_analysis(pseudo_op,
2699                                    xlim=(-2, 2), ylim=(-2, 2),
2700                                    xi_range=(0.1, 5), eta_range=(-5, 5),
2701                                    density=100):
2702        """
2703        Launch an interactive dashboard for symbol exploration using ipywidgets.
2704    
2705        This function provides a user-friendly interface to visualize various aspects of the pseudo-differential operator's symbol.
2706        It supports multiple visualization modes in both 1D and 2D, including group velocity fields, micro-support estimates,
2707        symplectic vector fields, symbol amplitude/phase, cotangent fiber structure, characteristic sets and Hamiltonian flows.
2708    
2709        Parameters
2710        ----------
2711        pseudo_op : PseudoDifferentialOperator
2712            The pseudo-differential operator whose symbol is to be analyzed interactively.
2713        xlim, ylim : tuple of float
2714            Spatial domain limits along x and y axes respectively.
2715        xi_range, eta_range : tuple
2716            Frequency domain limits along ξ and η axes respectively.
2717        density : int
2718            Number of points per axis used to construct the evaluation grid. Controls resolution.
2719    
2720        Notes
2721        -----
2722        - In 1D mode, sliders control the fixed frequency (ξ₀) and spatial position (x₀).
2723        - In 2D mode, additional sliders control the second frequency component (η₀) and second spatial coordinate (y₀).
2724        - Visualization updates dynamically as parameters are adjusted via sliders or dropdown menus.
2725        - Supported visualization modes:
2726            'Symbol Amplitude'           : |p(x,ξ)| or |p(x,y,ξ,η)|
2727            'Symbol Phase'               : arg(p(x,ξ)) or similar in 2D
2728            'Micro-Support (1/|p|)'      : Reciprocal of symbol magnitude
2729            'Cotangent Fiber'            : Structure of symbol over frequency space at fixed x
2730            'Characteristic Set'         : Zero set approximation {p ≈ 0}
2731            'Characteristic Gradient'    : |∇p(x, ξ)| or |∇p(x₀, y₀, ξ, η)|
2732            'Group Velocity Field'       : ∇_ξ p(x,ξ) or ∇_{ξ,η} p(x,y,ξ,η)
2733            'Symplectic Vector Field'    : (∇_ξ p, -∇_x p) or similar in 2D
2734            'Hamiltonian Flow'           : Trajectories generated by the Hamiltonian vector field
2735    
2736        Raises
2737        ------
2738        NotImplementedError
2739            If the spatial dimension is not 1D or 2D.
2740    
2741        Prints
2742        ------
2743        Interactive matplotlib figures with dynamic updates based on widget inputs.
2744        """
2745        dim = pseudo_op.dim
2746        expr = pseudo_op.expr
2747        vars_x = pseudo_op.vars_x
2748    
2749        mode_selector_1D = Dropdown(
2750            options=[
2751                'Symbol Amplitude',
2752                'Symbol Phase',
2753                'Micro-Support (1/|p|)',
2754                'Cotangent Fiber',
2755                'Characteristic Set',
2756                'Characteristic Gradient',
2757                'Group Velocity Field',
2758                'Symplectic Vector Field',
2759                'Hamiltonian Flow',
2760            ],
2761            value='Symbol Amplitude',
2762            description='Mode:'
2763        )
2764
2765        mode_selector_2D = Dropdown(
2766            options=[
2767                'Symbol Amplitude',
2768                'Symbol Phase',
2769                'Micro-Support (1/|p|)',
2770                'Cotangent Fiber',
2771                'Characteristic Set',
2772                'Characteristic Gradient',
2773                'Symplectic Vector Field',
2774                'Hamiltonian Flow',
2775            ],
2776            value='Symbol Amplitude',
2777            description='Mode:'
2778        )
2779    
2780        x_vals = np.linspace(*xlim, density)
2781        if dim == 2:
2782            y_vals = np.linspace(*ylim, density)
2783    
2784        if dim == 1:
2785            x, = vars_x
2786            xi = symbols('xi', real=True)
2787            grad_func = lambdify((x, xi), diff(expr, xi), 'numpy')
2788            symplectic_func = lambdify((x, xi), [diff(expr, xi), -diff(expr, x)], 'numpy')
2789            symbol_func = lambdify((x, xi), expr, 'numpy')
2790
2791            xi_slider = FloatSlider(min=xi_range[0], max=xi_range[1], step=0.1, value=1.0, description='ξ₀')
2792            x_slider = FloatSlider(min=xlim[0], max=xlim[1], step=0.1, value=0.0, description='x₀')
2793    
2794            def plot_1d(mode, xi0, x0):
2795                X = x_vals[:, None]
2796    
2797                if mode == 'Group Velocity Field':
2798                    V = grad_func(X, xi0)
2799                    plt.quiver(X, V, np.ones_like(V), V, scale=10, width=0.004)
2800                    plt.xlabel('x')
2801                    plt.title(f'Group Velocity Field at ξ={xi0:.2f}')
2802    
2803                elif mode == 'Micro-Support (1/|p|)':
2804                    Z = 1 / (np.abs(symbol_func(X, xi0)) + 1e-10)
2805                    plt.plot(x_vals, Z)
2806                    plt.xlabel('x')
2807                    plt.title(f'Micro-Support (1/|p|) at ξ={xi0:.2f}')
2808    
2809                elif mode == 'Symplectic Vector Field':
2810                    U, V = symplectic_func(X, xi0)
2811                    plt.quiver(X, V, U, V, scale=10, width=0.004)
2812                    plt.xlabel('x')
2813                    plt.title(f'Symplectic Field at ξ={xi0:.2f}')
2814    
2815                elif mode == 'Symbol Amplitude':
2816                    Z = np.abs(symbol_func(X, xi0))
2817                    plt.plot(x_vals, Z)
2818                    plt.xlabel('x')
2819                    plt.title(f'Symbol Amplitude |p(x,ξ)| at ξ={xi0:.2f}')
2820    
2821                elif mode == 'Symbol Phase':
2822                    Z = np.angle(symbol_func(X, xi0))
2823                    plt.plot(x_vals, Z)
2824                    plt.xlabel('x')
2825                    plt.title(f'Symbol Phase arg(p(x,ξ)) at ξ={xi0:.2f}')
2826    
2827                elif mode == 'Cotangent Fiber':
2828                    pseudo_op.visualize_fiber(x_vals, np.linspace(*xi_range, density), x0=x0)
2829    
2830                elif mode == 'Characteristic Set':
2831                    pseudo_op.visualize_characteristic_set(x_vals, np.linspace(*xi_range, density), x0=x0)
2832    
2833                elif mode == 'Characteristic Gradient':
2834                    pseudo_op.visualize_characteristic_gradient(x_vals, np.linspace(*xi_range, density), x0=x0)
2835    
2836                elif mode == 'Hamiltonian Flow':
2837                    pseudo_op.plot_hamiltonian_flow(x0=x0, xi0=xi0)
2838    
2839            # --- Dynamic container for sliders ---
2840            controls_box = VBox([mode_selector_1D, xi_slider, x_slider])
2841            # --- Function to adjust visible sliders based on mode ---
2842            def update_controls(change):
2843                mode = change['new']
2844                # modes that depend only on xi and eta
2845                if mode in ['Symbol Amplitude', 'Symbol Phase', 'Micro-Support (1/|p|)',
2846                            'Group Velocity Field', 'Symplectic Vector Field']:
2847                    controls_box.children = [mode_selector_1D, xi_slider]
2848                # modes that require xi and x
2849                elif mode in ['Hamiltonian Flow']:
2850                    controls_box.children = [mode_selector_1D, xi_slider, x_slider]
2851                # modes that require nothing
2852                elif mode in ['Cotangent Fiber', 'Characteristic Set', 'Characteristic Gradient']:
2853                    controls_box.children = [mode_selector_1D]
2854            mode_selector_1D.observe(update_controls, names='value')
2855            update_controls({'new': mode_selector_1D.value}) 
2856            # --- Interactive binding ---
2857            out = interactive_output(plot_1d, {'mode': mode_selector_1D, 'xi0': xi_slider, 'x0': x_slider})
2858            display(VBox([controls_box, out]))
2859
2860        elif dim == 2:
2861            x, y = vars_x
2862            xi, eta = symbols('xi eta', real=True)
2863            symplectic_func = lambdify((x, y, xi, eta), [diff(expr, xi), diff(expr, eta)], 'numpy')
2864            symbol_func = lambdify((x, y, xi, eta), expr, 'numpy')
2865
2866            xi_slider=FloatSlider(min=xi_range[0], max=xi_range[1], step=0.1, value=1.0, description='ξ₀')
2867            eta_slider=FloatSlider(min=eta_range[0], max=eta_range[1], step=0.1, value=1.0, description='η₀')
2868            x_slider=FloatSlider(min=xlim[0], max=xlim[1], step=0.1, value=0.0, description='x₀')
2869            y_slider=FloatSlider(min=ylim[0], max=ylim[1], step=0.1, value=0.0, description='y₀')
2870    
2871            def plot_2d(mode, xi0, eta0, x0, y0):
2872                X, Y = np.meshgrid(x_vals, y_vals, indexing='ij')
2873    
2874                if mode == 'Micro-Support (1/|p|)':
2875                    Z = 1 / (np.abs(symbol_func(X, Y, xi0, eta0)) + 1e-10)
2876                    plt.pcolormesh(X, Y, Z, shading='auto', cmap='inferno')
2877                    plt.colorbar(label='1/|p|')
2878                    plt.xlabel('x')
2879                    plt.ylabel('y')
2880                    plt.title(f'Micro-Support at ξ={xi0:.2f}, η={eta0:.2f}')
2881    
2882                elif mode == 'Symplectic Vector Field':
2883                    U, V = symplectic_func(X, Y, xi0, eta0)
2884                    plt.quiver(X, Y, U, V, scale=10, width=0.004)
2885                    plt.xlabel('x')
2886                    plt.ylabel('y')
2887                    plt.title(f'Symplectic Field at ξ={xi0:.2f}, η={eta0:.2f}')
2888    
2889                elif mode == 'Symbol Amplitude':
2890                    Z = np.abs(symbol_func(X, Y, xi0, eta0))
2891                    plt.pcolormesh(X, Y, Z, shading='auto')
2892                    plt.colorbar(label='|p(x,y,ξ,η)|')
2893                    plt.xlabel('x')
2894                    plt.ylabel('y')
2895                    plt.title(f'Symbol Amplitude at ξ={xi0:.2f}, η={eta0:.2f}')
2896    
2897                elif mode == 'Symbol Phase':
2898                    Z = np.angle(symbol_func(X, Y, xi0, eta0))
2899                    plt.pcolormesh(X, Y, Z, shading='auto', cmap='twilight')
2900                    plt.colorbar(label='arg(p)')
2901                    plt.xlabel('x')
2902                    plt.ylabel('y')
2903                    plt.title(f'Symbol Phase at ξ={xi0:.2f}, η={eta0:.2f}')
2904    
2905                elif mode == 'Cotangent Fiber':
2906                    pseudo_op.visualize_fiber(np.linspace(*xi_range, density), np.linspace(*eta_range, density),
2907                                              x0=x0, y0=y0)
2908    
2909                elif mode == 'Characteristic Set':
2910                    pseudo_op.visualize_characteristic_set(x_grid=x_vals, xi_grid=np.linspace(*xi_range, density),
2911                                                  y_grid=y_vals, eta_grid=np.linspace(*eta_range, density), x0=x0, y0=y0)
2912    
2913                elif mode == 'Characteristic Gradient':
2914                    pseudo_op.visualize_characteristic_gradient(x_grid=x_vals, xi_grid=np.linspace(*xi_range, density),
2915                                                  y_grid=y_vals, eta_grid=np.linspace(*eta_range, density), x0=x0, y0=y0)
2916    
2917                elif mode == 'Hamiltonian Flow':
2918                    pseudo_op.plot_hamiltonian_flow(x0=x0, y0=y0, xi0=xi0, eta0=eta0)
2919                    
2920            # --- Dynamic container for sliders ---
2921            controls_box = VBox([mode_selector_2D, xi_slider, eta_slider, x_slider, y_slider])
2922            # --- Function to adjust visible sliders based on mode ---
2923            def update_controls(change):
2924                mode = change['new']
2925                # modes that depend only on xi
2926                if mode in ['Symbol Amplitude', 'Symbol Phase', 'Micro-Support (1/|p|)', 'Symplectic Vector Field']:
2927                    controls_box.children = [mode_selector_2D, xi_slider, eta_slider]
2928                # modes that require xi, eta, x and y
2929                elif mode in ['Hamiltonian Flow']:
2930                    controls_box.children = [mode_selector_2D, xi_slider, eta_slider, x_slider, y_slider]
2931                # modes that require x and y
2932                elif mode in ['Cotangent Fiber', 'Characteristic Set', 'Characteristic Gradient']:
2933                    controls_box.children = [mode_selector_2D, x_slider, y_slider]
2934            mode_selector_2D.observe(update_controls, names='value')
2935            update_controls({'new': mode_selector_2D.value}) 
2936            # --- Interactive binding ---
2937            out = interactive_output(plot_2d, {'mode': mode_selector_2D, 'xi0': xi_slider, 'eta0': eta_slider, 'x0': x_slider, 'y0': y_slider})
2938            display(VBox([controls_box, out]))